Merge pull request #63 from SSCHAcode/new_cluster

A new cluster interface.

The new cluster interface waits for the queue without remaining pending and uses tar archives to move files between the server and the client, minimizing I/O operations.

Author: mesonepigreco

Modules/Cluster.py

@@ -1,6 +1,7 @@
 # -*- coding: utf-8 -*-
 from __future__ import print_function
 import sys, os
+import warnings
 import numpy as np
 import time
 
@@ -33,6 +34,9 @@
 
 import sscha.Parallel as Parallel
 from sscha.Parallel import pprint as print
+from sscha.Tools import NumpyEncoder
+
+import json
 
 import difflib
 
@@ -95,7 +99,11 @@
 UNITS_DEFAULT = "default"
 UNITS_HARTREE = "hartree"
 SUPPORTED_UNITS = [UNITS_DEFAULT, UNITS_HARTREE]
-
+class NumpyEncoder(json.JSONEncoder):
+    def default(self, obj):
+        if isinstance(obj, np.ndarray):
+            return obj.tolist()
+        return json.JSONEncoder.default(self, obj)
 
 class Ensemble:
     __debug_index__ = 0
@@ -179,7 +187,7 @@ def __init__(self, dyn0, T0, supercell = None, **kwargs):
         self.stress_computed = None
 
         # Get the extra quantities
-        self.extra_quantities = {}
+        self.all_properties = []
 
 
         # Setup the attribute control
@@ -610,7 +618,7 @@ def load_from_calculator_output(self, directory, out_ext = ".pwo"):
 
             # Get the stress if any
             try:
-                stress = ase_struct.get_stress(voigt=False)
+                stress = - ase_struct.get_stress(voigt=False) 
                 # eV/A^3 -> Ry/bohr^3
                 stress /= Rydberg / Bohr**3
                 count_stress += 1
@@ -762,6 +770,10 @@ def save_bin(self, data_dir, population_id = 1):
                 np.save("%s/stresses_pop%d.npy" % (data_dir, population_id), self.stresses)
 
             self.dyn_0.save_qe("%s/dyn_gen_pop%d_" % (data_dir, population_id))
+
+            if len(self.all_properties):
+                with open(os.path.join(data_dir, "all_properties_pop%d.json" % population_id), "w") as fp:
+                    json.dump({"properties" : self.all_properties}, fp, cls=NumpyEncoder)
 
     def save_enhanced_xyz(self, filename, append_mode = True, stress_key = "virial", forces_key = "force", energy_key = "energy"):
         """
@@ -861,7 +873,8 @@ def save_raw(self, root_directory, type_dict = None):
             np.savetxt(os.path.join(root_directory, "coord.raw"), self.xats.reshape((self.N, 3 * nat)))
 
             # Save the box
-            np.savetxt(os.path.join(root_directory, "box.raw"), np.tile(self.current_dyn.structure.unit_cell.ravel(), (self.N, 1)))
+            #Prepare an array with all the unit cells
+            np.savetxt(os.path.join(root_directory, "box.raw"), np.array([s.unit_cell.ravel() for s in self.structures]))
 
             # Save the forces
             np.savetxt(os.path.join(root_directory, "force.raw"), self.forces.reshape((self.N, 3*nat)) * Rydberg)
@@ -949,6 +962,16 @@ def load_bin(self, data_dir, population_id = 1, avoid_loading_dyn = False):
         # 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)
+
+        all_prop_fname = os.path.join(data_dir, "all_properties_pop%d.json" % population_id)
+        if os.path.exists(all_prop_fname):
+            with open(os.path.join(data_dir, "all_properties_pop%d.json" % population_id), "w") as fp:
+                props= json.load(fp)
+                if "properties" in props:
+                    self.all_properties = props["properties"]
+                else:
+                    warnings.warn("WARNING: found file {} but not able to load the properties keyword.".format(all_prop_fname))
+                    
 
 
     def init_from_structures(self, structures):
@@ -1616,6 +1639,39 @@ def get_free_energy(self, return_error = False):
             return free_energy, error
         return free_energy
 
+    def get_entropy(self, return_error = False):
+        r"""
+        GET THE ENTROPY
+        ===============
+
+        Get the total anharmonic entropy.
+
+        The equation implemented is the analytical derivative of the free energy,
+        and assumes that the SSCHA free energy is minimized.
+
+        .. math::
+
+            S = - \frac{dF}{dT}
+
+            S = S_{harm} - \left<V - {\mathcal V}\right>\sum_\mu \frac{1}{1 + 2n_\mu} \frac{dn_\mu}{dT}
+
+
+        where :math:`S_{harm}` is the 'harmonic' entropy computed from the dynamucal matrix, 
+        plus a correction accounting for the ensemble anharmonicity.
+
+        
+        Parameters
+        ----------
+            return_error : bool
+                If true, returns also the error
+
+        Results
+        -------
+            entropy, error : float
+                Returns the entropy and [optionally] the stochastic error.
+        """
+        raise NotImplementedError("Error, to be implemented")
+
 
     def get_free_energy_interpolating(self, target_supercell, support_dyn_coarse = None, support_dyn_fine = None, error_on_imaginary_frequency = True, return_error = False):
         """
@@ -2680,7 +2736,7 @@ def get_dynamical_bubble(self, q, w, smearing = 1e-5):
 
 
     def get_free_energy_hessian(self, include_v4 = False, get_full_hessian = True, verbose = False, \
-        use_symmetries = True):
+        use_symmetries = True, return_d3 = False):
         """
         GET THE FREE ENERGY ODD CORRECTION
         ==================================
@@ -2705,11 +2761,16 @@ def get_free_energy_hessian(self, include_v4 = False, get_full_hessian = True, v
             use_symmetries : bool
                 If true, the d3 and d4 are symmetrized in real space.
                 It requires that spglib is installed to detect symmetries in the supercell correctly.
+            return_d3 : bool
+                If true, returns also the tensor of three phonon scattering.
 
         Returns
         -------
             phi_sc : Phonons()
                 The dynamical matrix of the free energy hessian in (Ry/bohr^2)
+            d3 : ndarray (size = (3*nat_sc, 3*nat_sc, 3*nat_sc), Optional
+                Return the three-phonon-scattering tensor (in Ry atomic units).
+                Only if return_d3 is True. 
         """
         # For now the v4 is not implemented
         #     if include_v4:
@@ -2845,6 +2906,8 @@ def get_free_energy_hessian(self, include_v4 = False, get_full_hessian = True, v
                 dyn_hessian.dynmats[iq] = dynq_odd[iq, :, :] 
 
 
+        if return_d3:
+            return dyn_hessian, d3* 2.0 # Ha to Ry
         return dyn_hessian
 
 

Modules/Ensemble.py

@@ -46,6 +46,9 @@ def __init__(self, starting_unit_cell):
         self.uc_0_inv = np.linalg.inv(self.uc_0)
         self.x_start = None
         self.reset_strain = False
+        
+        # A control parameter for the line minimization (avoid failures)
+        self.last_factor = 0  
 
         self.uc_old = None
 
@@ -105,7 +108,8 @@ def get_line_step(self, grad):
             sys.stdout.flush()
 
             print("[CELL]  GRADIENT DOT DIRECTION = {}".format(factor))
-            
+            failure_lmin = False
+
             sys.stdout.flush()
             # Regularization (avoid run away)
             if factor > 2:
@@ -119,20 +123,36 @@ def get_line_step(self, grad):
                 print("                incrementing by a 15 %")
 
             elif factor < self.min_alpha_factor:
+                # Check if the last time was good
+                if self.last_factor != 0:
+                    if factor < self.last_factor:
+                        print("[CELL] Warning: reduced step and got worse; restarting line minimization with a good step.")
+                        failure_lmin = True
+                self.last_factor = factor
                 factor = self.min_alpha_factor
                 good_step = False
 
             if np.isnan(factor):
+                self.last_factor = 0
                 factor = self.min_alpha_factor
                 good_step = False
-                
-            self.alpha *= factor
+
+            # Check if the line minimization got contradditing result
+            # In this case stop the line minimization and proceed with a standard gradient descend.
+            if failure_lmin:
+                good_step = True
+            else:  
+                self.alpha *= factor
 
 
             #if self.alpha < self.min_alpha_step * self.alpha0:
             #    self.alpha = self.min_alpha_step * self.alpha0
             #self.alpha *= factor
 
+        # Reset the check on the line minimization if this is a good step
+        if good_step:
+            self.last_factor = 0
+
         return good_step
 
     def reset(self, unit_cell):
@@ -271,7 +291,7 @@ def UpdateCell(self, unit_cell, stress_tensor, fix_volume = False, verbose = Tru
 
 
 
-        # Get the strain tensor up to know
+        # Get the strain tensor up to now
         strain = np.transpose(self.uc_0_inv.dot(unit_cell) - I)
 
         # Get the gradient with respect to the strain        

Modules/Optimizer.py

@@ -441,6 +441,33 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
                 True if the minimization converged, False if the maximum number of 
                 populations has been reached.
         """
+
+        # Prepare the saving directory
+        if ensemble_loc is None:
+            ensemble_loc = self.data_dir
+
+        if (not ensemble_loc) and self.save_ensemble:
+            ERR_MSG = """
+Error, you must specify where to save the ensembles.
+       this can be done either passing ensemble_loc = "path/to/dir" 
+       for the ensemble, or by setting the data_dir attribute of this object.
+"""
+            raise IOError(ERR_MSG)
+        
+        if self.save_ensemble:
+            if not os.path.exists(ensemble_loc):
+                os.makedirs(ensemble_loc)
+            else:
+                if not os.path.isdir(ensemble_loc):
+                    ERR_MSG = """
+Error, the specified location to save the ensemble:
+       '{}' 
+       already exists and it is not a directory.
+""".format(ensemble_loc)
+                    raise IOError(ERR_MSG)
+
+        
+
         # Rescale the target pressure in eV / A^3
         target_press_evA3 = target_press / sscha.SchaMinimizer.__evA3_to_GPa__
         I = np.eye(3, dtype = np.float64)
@@ -509,6 +536,10 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
             self.minim.ensemble.dyn_0 = self.minim.dyn.Copy()
             if pop != start_pop or not restart_from_ens:
                 self.minim.ensemble.generate(self.N_configs)
+
+                # Save also the generation
+                #if ensemble_loc is not None and self.save_ensemble:
+                #    self.minim.ensemble.save_bin(ensemble_loc, pop)
 
                 # Compute energies and forces
                 self.minim.ensemble.compute_ensemble(self.calc, True, stress_numerical,

Modules/Relax.py

@@ -231,6 +231,7 @@ def __init__(self, ensemble = None, root_representation = "normal",
         self.__gw__ = []
         self.__gw_err__ = []
         self.__KL__ = []
+        self.__good_kasteps__ = []
 
 
 
@@ -554,6 +555,10 @@ def minimization_step(self, custom_function_gradient = None):
 
         # Update the previous gradient
         self.prev_grad = dyn_grad
+
+        # If the step is good, pass chose it
+        if self.minimizer.new_direction:
+            self.__good_kasteps__.append(len(self.__fe__) - 1)
 
     def setup_from_namelist(self, input_file):
         """
@@ -1232,7 +1237,25 @@ def finalize(self, verbose = 1):
                                                                     self.__gw_err__[-1] * __RyTomev__))
         print ("Kong-Liu effective sample size = ", self.ensemble.get_effective_sample_size())
         print ()
-        
+
+        # Print the total force on the structure
+        print("Total force on the centroids [eV/A]:")
+        forces, err = self.ensemble.get_average_forces(True)
+        # Apply the symmetries
+        if not self.neglect_symmetries:
+            qe_sym = CC.symmetries.QE_Symmetry(self.dyn.structure)
+            if self.use_spglib:
+                qe_sym.SetupFromSPGLIB()
+            else:
+                qe_sym.SetupQPoint()
+            qe_sym.SymmetrizeVector(forces)
+            qe_sym.SymmetrizeVector(err)
+            err /= np.sqrt(qe_sym.QE_nsym)
+
+        for i in range(self.dyn.structure.N_atoms):
+            print("{:4d}) {:14.6f}{:14.6f}{:14.6f}  +- {:14.6f}{:14.6f}{:14.6f}".format(*list([i] + list(forces[i, :] * CC.Units.RY_TO_EV) + list(err[i,:] * CC.Units.RY_TO_EV))))
+        print()
+
         if self.ensemble.has_stress and verbose >= 1:
             print ()
             print (" ==== STRESS TENSOR [GPa] ==== ")
@@ -1387,9 +1410,9 @@ def check_stop(self):
             print ()
             print ("The gw gradient satisfy the convergence condition.")
 
-        if self.gradi_op == "gc":
+        if self.gradi_op == "gc" or not self.minim_struct:
             total_cond = gc_cond
-        elif self.gradi_op == "gw":
+        elif self.gradi_op == "gw" or not self.minim_dyn:
             total_cond = gw_cond
         elif self.gradi_op == "all":
             total_cond = gc_cond and gw_cond
@@ -1465,18 +1488,18 @@ def plot_results(self, save_filename = None, plot = True):
 
             self.__gc__.append(self.__gc__[-1])
             self.__gc_err__.append(self.__gc_err__[-1])
+
+        fe = np.array(self.__fe__)[self.__good_kasteps__] * __RyTomev__
+        fe_err = np.array(self.__fe_err__)[self.__good_kasteps__] * __RyTomev__
 
-        # Convert the data in numpy arrays
-        fe = np.array(self.__fe__) * __RyTomev__
-        fe_err = np.array(self.__fe_err__) * __RyTomev__
 
-        gc = np.array(self.__gc__) * __RyTomev__
-        gc_err = np.array(self.__gc_err__) * __RyTomev__
+        gc = np.array(self.__gc__)[self.__good_kasteps__] * __RyTomev__
+        gc_err = np.array(self.__gc_err__)[self.__good_kasteps__] * __RyTomev__
 
-        gw = np.array(self.__gw__) * __RyTomev__
-        gw_err = np.array(self.__gw_err__) * __RyTomev__
+        gw = np.array(self.__gw__)[self.__good_kasteps__] * __RyTomev__
+        gw_err = np.array(self.__gw_err__)[self.__good_kasteps__] * __RyTomev__
 
-        kl = np.array(self.__KL__, dtype = np.float64)
+        kl = np.array(self.__KL__, dtype = np.float64)[self.__good_kasteps__]
 
         steps = np.arange(len(fe))
 
@@ -1489,6 +1512,7 @@ def plot_results(self, save_filename = None, plot = True):
             np.savetxt(save_filename, np.real(np.transpose(save_data)),
                        header = "Steps; Free energy +- error [meV]; FC gradient +- error [bohr^2]; Structure gradient +- error [meV / A]; Kong-Liu N_eff")
             print ("Minimization data saved in %s." % save_filename)
+            print('Good steps at {}'.format(self.__good_kasteps__))
 
 
         # Plot

Modules/SchaMinimizer.py

@@ -27,6 +27,12 @@
 import scipy, scipy.linalg
 import numpy as np
 import os
+import json
+class NumpyEncoder(json.JSONEncoder):
+    def default(self, obj):
+        if isinstance(obj, np.ndarray):
+            return obj.tolist()
+        return json.JSONEncoder.default(self, obj)
 
 
 # Here some usefull functions to solve linear systems