Merge branch 'master' into memopt-get-odd-straight

Author: mesonepigreco

.github/workflows/build.yml

@@ -20,7 +20,7 @@ jobs:
 
     services:
       rabbitmq:
-          image: rabbitmq:latest
+          image: rabbitmq:3.12-management
           ports:
           -   5672:5672
 

.github/workflows/python-testsuite.yml

@@ -333,6 +333,30 @@ def __setattr__(self, name, value):
             super(Cluster, self).__setattr__(name, value)
 
 
+    def __getstate__(self):
+        """
+        Return the picklable state of the cluster.
+
+        The thread lock created by compute_ensemble_batch cannot be pickled,
+        so it is dropped here. This allows sscha.Utilities.save_binary to
+        store objects holding a cluster after a calculation has run.
+        """
+        state = self.__dict__.copy()
+        state["lock"] = None
+        return state
+
+
+    def __setstate__(self, state):
+        """
+        Restore the cluster from a pickled state.
+
+        The thread lock is transient runtime state and is reset to None,
+        as after __init__; compute_ensemble_batch recreates it when needed.
+        """
+        state["lock"] = None
+        self.__dict__.update(state)
+
+
 
     def copy_file(self, source, destination, server_source = False, server_dest = True, raise_error=False, **kwargs):
         """

Modules/Cluster.py

@@ -4230,7 +4230,57 @@ def w_to_a(self,w, T):
             a[:] = np.sqrt((1.0 / np.tanh(0.5 * w * 315774.65221921849 / T)) / (2.0 * w))
         return a
 
+
+
+
+# ------------------------------------------------------------------------------
+def load_ensemble_bin(directory, population, temperature, nqirr=-1):
+    """
+    Load the ensemble from the given directory. This subroutine automatically initializes the ensemble
+    and avoids the hursle of having to define and initialize a dynamical matrix before loading the ensemble.
+
+
+    Parameters
+    ----------
+        - directory : string
+            Path to the directory where the ensemble is located
+        - population : Int
+            The index of the ensemble
+        - temperature : Float
+            The temperature for the ensemble
+        - nqirr : Int
+            The number of irreducible points of the dynamical matrix.
+            If not provided (or negative), it is inferred from the files.
+
+    Returns
+    -------
+        - ensemble : Ensemble()
+            The Ensemble object.
+    """
+
+    # Get the generated dyn
+    generated_name = os.path.join(directory, "dyn_gen_pop%d_" % population)
+
+    # Infer the nqirr if not provided
+    if nqirr < 0:
+        nqirr = len([x for x in os.listdir(directory) if x.startswith("dyn_gen_pop%d_" % population)])
+        if os.path.exists(generated_name + "0"):
+            nqirr -= 1
+
+    # Load the dynamical matrix
+    dyn = CC.Phonons.Phonons(generated_name, nqirr)
+
+    # Check if nqirr was correct
+    assert np.prod(dyn.GetSupercell()) == len(dyn.dynmats), "Error, the value of nqirr = {} is wrong. Check wether you specified the correct value.".format(nqirr)
+
+    ensemble = Ensemble(dyn, temperature)
+    ensemble.load_bin(directory, population)
+
+    return ensemble
+
+
 #-------------------------------------------------------------------------------
+
 def _wrapper_julia_get_upsilon_q(*args, **kwargs):
     """Worker function, just for testing"""
     return julia.Main.get_upsilon_fourier(*args, **kwargs)

Modules/Ensemble.py

@@ -224,6 +224,9 @@ def run_step(self, gradient, kl_new):
             # Enlarge the step
             if not self.fixed_step:
                 self.step *= self.increment_step
+            
+            # Perform the minimization step for new direction
+            self.current_x = self.old_x - self.step * self.direction
         else:
             # Proceed with the line minimization
 
@@ -243,14 +246,16 @@ def run_step(self, gradient, kl_new):
                     print("Step too large (scalar = {} | kl_ratio = {}), reducing to {}".format(scalar, kl_ratio, self.step))
                     #print("Direction: ", self.direction)
                     #print("Gradient: ", gradient)
+                
+                # Try again with reduced step
+                self.current_x = self.old_x - self.step * self.direction
             else:
                 # The step is good, therefore next step perform a new direction
                 self.new_direction = True
                 if self.verbose:
                     print("Good step found with {}, try increment".format(self.step))
-
-        # Perform the minimiziation step
-        self.current_x = self.old_x - self.step * self.direction
+                # DO NOT update current_x - we accept the current position
+                # (current_x was already updated in the previous step)
 
 
     def update_dyn(self, new_kl_ratio, dyn_gradient, structure_gradient = None):

Modules/Minimizer.py

@@ -634,7 +634,7 @@ def minimization_step(self, custom_function_gradient = None, timer=None):
 
             # Update the structure
             if self.minim_struct:
-                self.dyn.structure.coords[:,:] = new_struct
+                self.dyn.structure.coords[:,:] = np.real(new_struct)
 
             # Check if we must enforce the symmetries and the sum rule:
             if self.enforce_sum_rule and (not self.neglect_symmetries):
@@ -643,6 +643,11 @@ def minimization_step(self, custom_function_gradient = None, timer=None):
                 else:
                     self.dyn.Symmetrize(use_spglib = self.use_spglib)
 
+                # Enforce the dynamical matrix to be real at q = -q + G (time-reversal symmetry)
+                bg = self.dyn.structure.get_reciprocal_vectors() / (2 * np.pi)
+                for iq, q in enumerate(self.dyn.q_tot):
+                    if CC.Methods.get_min_dist_into_cell(bg, q, -q) < 1e-6:
+                        self.dyn.dynmats[iq] = np.real(self.dyn.dynmats[iq])
 
             # If we have imaginary frequencies, force the kl ratio to zero
 
@@ -964,11 +969,11 @@ def print_info(self):
         print (" supercell size = ", " ".join([str(x) for x in self.ensemble.supercell]))
 
         # Get the current frequencies
-        w, pols = self.dyn.GenerateSupercellDyn(self.ensemble.supercell).DyagDinQ(0)
+        w, pols = self.dyn.DiagonalizeSupercell()#self.dyn.GenerateSupercellDyn(self.ensemble.supercell).DyagDinQ(0)
         w *= __RyToCm__
 
         # Get the starting frequencies
-        w0, p0 = self.ensemble.dyn_0.GenerateSupercellDyn(self.ensemble.supercell).DyagDinQ(0)
+        w0, p0 = self.ensemble.dyn_0.DiagonalizeSupercell()
         w0 *= __RyToCm__
 
         print ()

Modules/SchaMinimizer.py

@@ -5,4 +5,4 @@
 """
 
 
-__all__ = ["Ensemble", "SchaMinimizer"]
+__all__ = ["Ensemble", "SchaMinimizer", "cli"]