Merge branch 'master' into examples/aiida

Author: mesonepigreco

Modules/Relax.py

@@ -584,8 +584,6 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
             if kind_minimizer == "RPSD":
                 # Compute the static bulk modulus
                 sbm = GetStaticBulkModulus(self.minim.dyn.structure, self.calc)
-                print ("BM:")
-                print (sbm)
                 BFGS = sscha.Optimizer.SD_PREC_UC(self.minim.dyn.structure.unit_cell, sbm)
 
             # Generate the ensemble
@@ -607,16 +605,13 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
                         self.calc, True, stress_numerical, cluster = self.cluster)
                 #self.minim.ensemble.get_energy_forces(self.calc, True, stress_numerical = stress_numerical)
 
-                print("RELAX force length:", len(self.minim.ensemble.force_computed))
 
                 if ensemble_loc is not None and self.save_ensemble:
                     self.minim.ensemble.save_bin(ensemble_loc, pop)
-                print("RELAX force length:", len(self.minim.ensemble.force_computed))
 
             self.minim.population = pop
             self.minim.init(delete_previous_data = False)
 
-            print("RELAX force length:", len(self.minim.ensemble.force_computed))
             self.minim.run(custom_function_pre = self.__cfpre__,
                            custom_function_post = self.__cfpost__,
                            custom_function_gradient = self.__cfg__)
@@ -633,6 +628,14 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
             # Get the pressure
             Press = np.trace(stress_tensor) / 3
 
+            # Get the error correctly accounting for symmetries
+            stress_diagonal_err = np.diag(stress_err)
+            press_err = 0
+            for i in range(3):
+                if not stress_diagonal_err[i] in stress_diagonal_err[:i]:
+                    press_err += stress_diagonal_err[i]**2
+            press_err = np.sqrt(press_err)
+
             # Get the volume
             Vol = self.minim.dyn.structure.get_volume()
 
@@ -660,6 +663,11 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
  ENTHALPIC CONTRIBUTION
  ======================
 
+ Current pressure P = {:.4f} +- {:.4f} GPa
+ Target pressure  P = {:.4f} GPa
+
+ For enthalpy we use the target pressure.
+
  P = {:.4f} GPa   V = {:.4f} A^3
 
  P V = {:.8e} eV
@@ -668,7 +676,7 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
  Gibbs Free energy = {:.10e} eV {}
  Zero energy = {:.10e} eV
 
- """.format(target_press , Vol,target_press_evA3 * Vol, helmoltz, mark_helmoltz, gibbs, mark_gibbs, self.minim.eq_energy)
+ """.format(Press / CC.Units.GPA_TO_EV_PER_A3, press_err/ CC.Units.GPA_TO_EV_PER_A3, target_press , Vol,target_press_evA3 * Vol, helmoltz, mark_helmoltz, gibbs, mark_gibbs, self.minim.eq_energy)
             print(message)
             # print " ====================== "
             # print " ENTHALPIC CONTRIBUTION "