Merge pull request #185 from ICAMS/get_energy_dissipation

Get energy dissipation

Author: srmnitc

calphy/integrators.py

@@ -226,6 +226,7 @@ def integrate_rs(simfolder, f0, t,
     
     """
     ws = []
+    es = []
     p = p/(10000*160.21766208)
 
     for i in range(1, nsims+1):
@@ -245,9 +246,12 @@ def integrate_rs(simfolder, f0, t,
         wf = cumtrapz(fdx, flambda,initial=0)
         wb = cumtrapz(bdx[::-1], blambda[::-1],initial=0)
         w = (wf + wb) / (2*flambda)
+        e = np.max(np.abs((wf - wb)/(2*flambda)))
+
         ws.append(w)
-    
+        es.append(e)
 
+    e_diss = np.min(es)
     wmean = np.mean(ws, axis=0)
     werr = np.std(ws, axis=0)
     temp = t/flambda
@@ -257,8 +261,9 @@ def integrate_rs(simfolder, f0, t,
     if not return_values:
         outfile = os.path.join(simfolder, "temperature_sweep.dat")
         np.savetxt(outfile, np.column_stack((temp, f, werr)))
+        return None, e_diss
     else:
-        return (temp, f, werr)
+        return (temp, f, werr), e_diss
 
 
 def integrate_ps(simfolder, f0, natoms, pi, pf, nsims=1, 

calphy/phase.py

@@ -1206,7 +1206,7 @@ def integrate_reversible_scaling(self, scale_energy=True, return_values=False):
         res : list of lists of shape 1x3
             Only returned if `return_values` is True.
         """
-        res = integrate_rs(
+        res, ediss = integrate_rs(
             self.simfolder,
             self.fe,
             self.calc._temperature,
@@ -1217,6 +1217,11 @@ def integrate_reversible_scaling(self, scale_energy=True, return_values=False):
             return_values=return_values,
         )
 
+        self.logger.info(f'Maximum energy dissipation along the temperature scaling part: {ediss} eV/atom')
+        if np.abs(ediss) > 1E-4:
+            self.logger.warning(f'Found max energy dissipation of {ediss} along the temperature scaling path. Please ensure there are no structural changes!')
+
+
         if return_values:
             return res