Merge branch 'line_minimization' of github.com:SSCHAcode/python-sscha into line_minimization

Author: mesonepigreco

Modules/Ensemble.py

@@ -836,13 +836,11 @@ def load_bin(self, data_dir, population_id = 1, avoid_loading_dyn = False):
 
         # Setup the initial weights
         self.rho = np.ones(self.N, dtype = np.float64)
-
-
+        
         # Setup that both forces and stresses are not computed
         self.stress_computed = np.ones(self.N, dtype = bool)
         self.force_computed = np.ones(self.N, dtype = bool)
 
-        
 
     def init_from_structures(self, structures):
         """

Modules/Relax.py

@@ -98,6 +98,12 @@ def __init__(self, minimizer = None, ase_calculator=None, N_configs=1, max_pop =
         self.target_pressure = 0
         self.fix_volume = False
 
+        # Options for the cell relaxation
+        # If true the cell shape is fixed in a variable-cell relaxation
+        # even if the symmetries allows for more degrees of freedom in the cell shape.
+        # Usefull (for example) if you want to enfoce a cubic cell even if the structure brakes the symmetries
+        self.fix_cell_shape = False  
+
 
         # Setup the attribute control
         self.__total_attributes__ = [item for item in self.__dict__.keys()]
@@ -351,7 +357,10 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
         This function performs a variable cell SCHA relaxation at constant pressure,
         It is similar to the relax calculation, but the unit cell is updated according
         to the anharmonic stress tensor at each new population. 
-        The cell optimization is performed with the BFGS algorithm. 
+
+        By default, all the degrees of freedom compatible with the symmetry group are relaxed in the cell.
+        You can constrain the cell to keep the same shape by setting fix_cell_shape = True.
+
         
         NOTE: 
             remember to setup the stress_offset variable of the SCHA_Minimizer,
@@ -360,6 +369,7 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
             stress tensor difference between a single very high-cutoff calculation and a
             single low-cutoff (the one you use), this difference will be added at the final
             stress tensor to get a better estimation of the true stress.
+
         
         Parameters
         ----------
@@ -393,7 +403,8 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
                 does not support it by default)
             cell_relax_algorithm : string
                 This identifies the stress algorithm. It can be both sd (steepest-descent),
-                cg (conjugate-gradient) or bfgs (Quasi-newton)
+                cg (conjugate-gradient) or bfgs (Quasi-newton).
+                The most robust one is SD. Do not change if you are not sure what you are doing.
             fix_volume : bool, optional
                 If true (default False) the volume is fixed, therefore only the cell shape is relaxed.
                 
@@ -558,7 +569,11 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
             # print ""
 
             # Perform the cell step
-            cell_gradient = (stress_tensor - I *target_press_evA3)
+            if self.fix_cell_shape:
+                # Use a isotropic stress tensor to keep the same cell shape
+                cell_gradient = I * (Press - target_press_evA3)
+            else:
+                cell_gradient = (stress_tensor - I *target_press_evA3)
 
             new_uc = self.minim.dyn.structure.unit_cell.copy()
             BFGS.UpdateCell(new_uc,  cell_gradient, fix_volume)

Modules/SchaMinimizer.py

@@ -160,6 +160,8 @@ def __init__(self, ensemble = None, root_representation = "normal",
 
         self.minimizer = None# 
 
+        self.max_diag_error_counter = __MAX_DIAG_ERROR_COUNTER__
+
         # Projection. This is chosen to fix some constraint on the minimization
         self.projector_dyn = None
         self.projector_struct = None
@@ -190,6 +192,9 @@ def __init__(self, ensemble = None, root_representation = "normal",
         # If true use spglib to impose symmetries
         # NOTE: It is slower, but it enables using supercells
         self.use_spglib = False
+
+        # If True, enforce the symmetrization and the sum rule after each step
+        self.enforce_sum_rule = False
 
 
         # Setup the statistical threshold
@@ -461,32 +466,41 @@ def minimization_step(self, custom_function_gradient = None):
             if self.minim_struct:
                 self.dyn.structure.coords[:,:] = new_struct
 
+            # Check if we must enforce the symmetries and the sum rule:
+            if self.enforce_sum_rule and (not self.neglect_symmetries):
+                self.dyn.Symmetrize(use_spglib = self.use_spglib)
+
+
             # If we have imaginary frequencies, force the kl ratio to zero
             if self.check_imaginary_frequencies():
                 print("Immaginary frequencies found! Redoing the step.")
                 new_kl_ratio = 0
                 is_diag_ok = False
-            
-
-            # Update the ensemble
-            try:
-                self.update()
-            except np.linalg.LinAlgError as error:
-                print("Diagonalization error:")
-                print(error)
-                print("Reducing the minimization step...")
-                new_kl_ratio = 0 # Force step reduction
-                is_diag_ok = False 
-                diag_error_counter += 1
-            
-            if diag_error_counter >= __MAX_DIAG_ERROR_COUNTER__:
+                diag_error_counter += 1            
+            else:
+                # Update the ensemble
+                try:
+                    self.update()
+                except np.linalg.LinAlgError as error:
+                    print("Diagonalization error:")
+                    print(error)
+                    print("Reducing the minimization step...")
+                    new_kl_ratio = 0 # Force step reduction
+                    is_diag_ok = False 
+                    diag_error_counter += 1
+            
+            if diag_error_counter >= self.max_diag_error_counter:
                 ERROR_MSG = """
-Error, exceeded the maximum number of diagonalization error. 
-       something is very wrong with the dynamical matrix.
+Error, exceeded the maximum number of diagonalization error. (or imaginary steps)
+
+       you can get rid of this error by increasing max_diag_error_counter variable.
+
+       Something is very wrong with the dynamical matrix.
        I'm saving the dynamical matrix as error_dyn, 
        if you want to check it or restart the calculation from there.
 """
                 print(ERROR_MSG)
+                sys.stdout.flush()
                 self.dyn.save_qe("error_dyn")
                 raise ValueError(ERROR_MSG)
 
@@ -1227,11 +1241,14 @@ def check_imaginary_frequencies(self):
         """
 
         # Get the frequencies
-        superdyn = self.dyn.GenerateSupercellDyn(self.ensemble.supercell)
-        w, pols = superdyn.DyagDinQ(0)
+        #superdyn = self.dyn.GenerateSupercellDyn(self.ensemble.supercell)
+        w, pols = self.dyn.DiagonalizeSupercell()#.DyagDinQ(0)
+        ss = self.dyn.structure.generate_supercell(self.dyn.GetSupercell())
 
         # Get translations
-        trans_mask = ~CC.Methods.get_translations(pols, superdyn.structure.get_masses_array())
+        trans_mask = ~CC.Methods.get_translations(pols, ss.get_masses_array())
+
+        current_n_trans = np.sum((~trans_mask).astype(int))
 
         # Remove translations
         w = w[trans_mask]
@@ -1240,10 +1257,30 @@ def check_imaginary_frequencies(self):
         # Frequencies are ordered, check if the first one is negative.
         if w[0] < 0:
             print ("Total frequencies (excluding translations):")
-            superdyn0 = self.ensemble.dyn_0.GenerateSupercellDyn(self.ensemble.supercell)
-            wold, pold = superdyn0.DyagDinQ(0)
+            #superdyn0 = self.ensemble.dyn_0.GenerateSupercellDyn(self.ensemble.supercell)
+            wold, pold = self.ensemble.dyn_0.DiagonalizeSupercell()# superdyn0.DyagDinQ(0)
+
+            ss0 = self.ensemble.dyn_0.structure.generate_supercell(self.dyn.GetSupercell())
 
-            trans_mask = ~CC.Methods.get_translations(pold, superdyn0.structure.get_masses_array())
+            trans_mask = ~CC.Methods.get_translations(pold, ss0.get_masses_array())
+
+            old_n_trans = np.sum((~trans_mask).astype(int))
+
+            if old_n_trans != current_n_trans:
+                ERR_MSG = """
+Error, the number of translational modes before and after the step is different.
+    original translational modes: {}
+    new translational modes: {}
+
+You can try to fix this error setting the {} variable of {} class to True.
+    """.format(old_n_trans, current_n_trans, "enforce_sum_rule", self.__class__.__name__)
+                print(ERR_MSG)
+                if not self.enforce_sum_rule:
+                    raise ValueError(ERR_MSG)
+                else:
+                    return True
+
+
             wold = wold[trans_mask] * __RyToCm__
             pold = pold[:, trans_mask]
             total_mask = list(range(len(w)))