diff --git a/src/PERI/pair_peri_eps.cpp b/src/PERI/pair_peri_eps.cpp
index fe03ac53a1..0bf339933f 100644
--- a/src/PERI/pair_peri_eps.cpp
+++ b/src/PERI/pair_peri_eps.cpp
@@ -231,9 +231,8 @@ void PairPeriEPS::compute(int eflag, int vflag)
    
 
   // compute the dilatation on each particle
-
   compute_dilatation();
-
+  
   // communicate dilatation (theta) of each particle
   comm->forward_comm_pair(this);
 
@@ -347,7 +346,6 @@ void PairPeriEPS::compute(int eflag, int vflag)
               
       double deviatoric_extension = dr - (theta[i]* r0[i][jj] / 3.0);
       edpNp1 = deviatorPlasticextension[i][jj];
-
   
       double tdtrialValue = ( 15 * shearmodulus[itype][itype]) *
         ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * 
@@ -378,9 +376,7 @@ void PairPeriEPS::compute(int eflag, int vflag)
 
       // since I-J is double counted, set newton off & use 1/2 factor and I,I
 
-      if (eflag) evdwl =  (0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
-                    omega_plus * deviatoric_extension * deviatoric_extension) + 
-                    (0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
+      if (eflag) evdwl =  (0.5 * 15 * shearmodulus[itype][itype]/wvolume[i] *
                        omega_plus * (deviatoric_extension - edpNp1) * 
                       (deviatoric_extension-edpNp1)) * vfrac[j] * vfrac_scale;
       if (evflag) ev_tally(i,i,nlocal,0,0.5*evdwl,0.0,
@@ -402,7 +398,7 @@ void PairPeriEPS::compute(int eflag, int vflag)
                          (alpha[itype][jtype] * stretch));
 
       first = false;
-    }
+    }    
   }
 
   // store new s0
@@ -415,8 +411,7 @@ void PairPeriEPS::compute(int eflag, int vflag)
          double temp_data = deviatorPlasticExtTemp[i][jj];
          deviatorPlasticextension[i][jj] = temp_data;
       }
-  }    
-         
+  }            
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/PERI/pair_peri_eps.h b/src/PERI/pair_peri_eps.h
index c6eacce1c6..bf20becc5d 100644
--- a/src/PERI/pair_peri_eps.h
+++ b/src/PERI/pair_peri_eps.h
@@ -27,6 +27,7 @@ namespace LAMMPS_NS {
 class PairPeriEPS : public Pair {
  public:
   double *theta;
+  double *elastic_energy;
 
   PairPeriEPS(class LAMMPS *);
   virtual ~PairPeriEPS();
diff --git a/src/PERI/pair_peri_lps.cpp b/src/PERI/pair_peri_lps.cpp
index c8ee9548a5..cb17d026cf 100644
--- a/src/PERI/pair_peri_lps.cpp
+++ b/src/PERI/pair_peri_lps.cpp
@@ -285,6 +285,7 @@ void PairPeriLPS::compute(int eflag, int vflag)
 
       omega_plus  = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0);
       omega_minus = influence_function(delx0,dely0,delz0);
+      
       rk = ( (3.0 * bulkmodulus[itype][itype]) -
              (5.0 * shearmodulus[itype][itype]) ) * vfrac[j] * vfrac_scale *
         ( (omega_plus * theta[i] / wvolume[i]) +
@@ -302,6 +303,8 @@ void PairPeriLPS::compute(int eflag, int vflag)
       // since I-J is double counted, set newton off & use 1/2 factor and I,I
 
       double deviatoric_extension = dr - (theta[i]* r0[i][jj] / 3.0);
+                  
+      
       if (eflag) evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
                    omega_plus*(deviatoric_extension * deviatoric_extension) *
                    vfrac[j] * vfrac_scale;
@@ -323,7 +326,7 @@ void PairPeriLPS::compute(int eflag, int vflag)
                          (alpha[itype][jtype] * stretch));
 
       first = false;
-    }
+    }      
   }
 
   // store new s0
@@ -618,6 +621,7 @@ void PairPeriLPS::compute_dilatation()
   }
 }
 
+
 /* ----------------------------------------------------------------------
    communication routines
  ---------------------------------------------------------------------- */
diff --git a/src/PERI/pair_peri_lps.h b/src/PERI/pair_peri_lps.h
index c3593e57ac..d28a2385b7 100644
--- a/src/PERI/pair_peri_lps.h
+++ b/src/PERI/pair_peri_lps.h
@@ -27,6 +27,7 @@ namespace LAMMPS_NS {
 class PairPeriLPS : public Pair {
  public:
   double *theta;
+  double *elastic_energy;
 
   PairPeriLPS(class LAMMPS *);
   virtual ~PairPeriLPS();
diff --git a/src/PERI/pair_peri_ves.cpp b/src/PERI/pair_peri_ves.cpp
index baa826ea82..6668de9a2e 100644
--- a/src/PERI/pair_peri_ves.cpp
+++ b/src/PERI/pair_peri_ves.cpp
@@ -217,12 +217,10 @@ void PairPeriVES::compute(int eflag, int vflag)
     memory->create(theta,nmax,"pair:theta");
   }
 
-  // compute the dilatation on each particle
-
+  // Compute the dilatation on each particle
   compute_dilatation();
 
   // communicate dilatation (theta) of each particle
-
   comm->forward_comm_pair(this);
 
   // communicate weighted volume (wvolume) upon every reneighbor
@@ -310,7 +308,6 @@ void PairPeriVES::compute(int eflag, int vflag)
       else fbond = 0.0;
 
       // for viscoelasticity
-
       lambdai=m_lambdai[itype][itype];
       double taui = m_taubi[itype][itype];  
       double c1 = taui/timestepsize;
@@ -349,12 +346,12 @@ void PairPeriVES::compute(int eflag, int vflag)
 
       // since I-J is double counted, set newton off & use 1/2 factor and I,I
 
-      if (eflag) evdwl =  (0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
+      if (eflag) evdwl =  ((0.5 * 15 * (1 - lambdai) * shearmodulus[itype][itype]/wvolume[i] *
                     omega_plus * deviatoric_extension * 
                     deviatoric_extension) + 
-                    (0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) *
+                    (0.5 * 15 * lambdai * shearmodulus[itype][itype]/wvolume[i] *
                     omega_plus * (deviatoric_extension-edbNp1) * 
-                    (deviatoric_extension-edbNp1)) * vfrac[j] * vfrac_scale;
+                    (deviatoric_extension-edbNp1))) * vfrac[j] * vfrac_scale;
       if (evflag) ev_tally(i,i,nlocal,0,0.5*evdwl,0.0,
                            0.5*fbond*vfrac[i],delx,dely,delz);
 
@@ -378,7 +375,7 @@ void PairPeriVES::compute(int eflag, int vflag)
                          (alpha[itype][jtype] * stretch));
 
       first = false;
-    }
+    }  
   }
 
   // store new s0
@@ -620,6 +617,7 @@ void PairPeriVES::compute_dilatation()
   double lc = domain->lattice->xlattice;
   double half_lc = 0.5*lc;
 
+
   double **r0   = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0;
   tagint **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner;
   int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner;
@@ -693,6 +691,7 @@ void PairPeriVES::compute_dilatation()
   }
 }
 
+
 /* ----------------------------------------------------------------------
    communication routines
 ---------------------------------------------------------------------- */
diff --git a/src/PERI/pair_peri_ves.h b/src/PERI/pair_peri_ves.h
index fbf6ab8bd9..4e3db61d0a 100644
--- a/src/PERI/pair_peri_ves.h
+++ b/src/PERI/pair_peri_ves.h
@@ -27,6 +27,7 @@ namespace LAMMPS_NS {
 class PairPeriVES : public Pair {
  public:
   double *theta;
+  double *elastic_energy;  
 
   PairPeriVES(class LAMMPS *);
   virtual ~PairPeriVES();
@@ -44,7 +45,8 @@ class PairPeriVES : public Pair {
   void read_restart_settings(FILE *) {}
   double memory_usage();
   double influence_function(double, double, double);
-  void compute_dilatation();
+  void compute_dilatation();  
+  
 
  protected:
   int ifix_peri;