Merge pull request #166 from SSCHAcode/warning_on_stress_vcrelax

Added an error if no stress tensor is found

Author: mesonepigreco

Modules/Relax.py

@@ -609,6 +609,8 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
                 if ensemble_loc is not None and self.save_ensemble:
                     self.minim.ensemble.save_bin(ensemble_loc, pop)
 
+
+
             self.minim.population = pop
             self.minim.init(delete_previous_data = False)
 
@@ -625,6 +627,11 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
             stress_tensor *= sscha.SchaMinimizer.__RyBohr3_to_evA3__
             stress_err *=  sscha.SchaMinimizer.__RyBohr3_to_evA3__
 
+            # Check if the stress tensor is actually loaded in the ensemble
+            if np.max(np.abs(self.minim.ensemble.stresses)) < 1e-10:
+                # Probably there is an error in the submission of the stress tensor calculation
+                raise ValueError("Error, the stress tensor is not loaded in the ensemble. Check the stress tensor calculation.")
+
             # Get the pressure
             Press = np.trace(stress_tensor) / 3
 
@@ -676,7 +683,7 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
  Gibbs Free energy = {:.10e} eV {}
  Zero energy = {:.10e} eV
 
- """.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)
+ """.format(Press / CC.Units.GPA_TO_EV_PER_A3, press_err/ CC.Units.GPA_TO_EV_PER_A3, target_press, target_press, Vol, target_press_evA3 * Vol, helmoltz, mark_helmoltz, gibbs, mark_gibbs, self.minim.eq_energy)
             print(message)
             # print " ====================== "
             # print " ENTHALPIC CONTRIBUTION "

tests/test_vc_relax/test_vcrelax_dummy.py

@@ -0,0 +1,65 @@
+import cellconstructor as CC, cellconstructor.Phonons
+import sscha, sscha.Ensemble
+import sscha.SchaMinimizer, sscha.Relax
+import sscha.Utilities
+
+import ase, ase.calculators, ase.calculators.emt
+import numpy as np
+import os, sys
+
+import pytest
+
+
+def test_gold(use_julia=True, verbose=False):
+    supercell=(2,2,2)
+    
+    # Change directory to the current script
+    total_path = os.path.dirname(os.path.abspath(__file__))
+    os.chdir(total_path)
+    
+    cif_file="../test_release/Au_mp-81_primitive.cif"
+    temperature = 300
+
+    # Load the structure
+    struct = CC.Structure.Structure()
+    struct.read_generic_file(cif_file)
+    gold_dyn = CC.Phonons.compute_phonons_finite_displacements(struct, 
+        ase.calculators.emt.EMT(), supercell=supercell)
+
+    # Setup a sscha
+    gold_dyn.Symmetrize()
+    gold_dyn.ForcePositiveDefinite()
+
+    ensemble = sscha.Ensemble.Ensemble(gold_dyn, temperature)
+    minim = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
+    minim.set_minimization_step(0.001)
+    minim.meaningful_factor = 1e-5
+    minim.use_julia = use_julia
+
+    relax = sscha.Relax.SSCHA(minim,
+            ase_calculator=ase.calculators.emt.EMT(),
+            N_configs = 50,
+            max_pop = 3)
+    relax.save_ensemble=False
+
+    cfp=None
+    if verbose:
+        ioinfo = sscha.Utilities.IOInfo()
+        ioinfo.SetupSaving("minim")
+        cfp = ioinfo.CFP_SaveAll
+
+    relax.setup_custom_functions(
+            custom_function_post=cfp)
+
+    relax.vc_relax(target_press=0)
+    
+    if verbose:
+        relax.minim.dyn.save_qe("final_dyn")
+
+if __name__ == "__main__":
+    use_julia=True
+    if len(sys.argv) > 1:
+        if sys.argv[1].lower() != "julia":
+            use_julia=False
+                
+    test_gold(verbose=True, use_julia=use_julia)