Merge pull request #83 from SebastianHavens/lammps-seed-bug

Ensure that random seed for initialising velocities is never zero
ICAMS · Sep 14, 2023 · 332e4bd · 332e4bd
2 parents f9aa07a + 8929b6c
commit 332e4bd
Show file tree

Hide file tree

Showing 4 changed files with 18 additions and 17 deletions.
diff --git a/calphy/alchemy.py b/calphy/alchemy.py
@@ -158,7 +158,7 @@ def run_integration(self, iteration=1):
         #remap the box to get the correct pressure
         lmp = ph.remap_box(lmp, self.lx, self.ly, self.lz)
 
-        lmp.command("velocity          all create %f %d mom yes rot yes dist gaussian"%(self.calc._temperature, np.random.randint(0, 10000)))
+        lmp.command("velocity          all create %f %d mom yes rot yes dist gaussian"%(self.calc._temperature, np.random.randint(1, 10000)))
         # Integrator & thermostat.
         if self.calc._npt:
             lmp.command("fix             f1 all npt temp %f %f %f %s %f %f %f"%(self.calc._temperature, self.calc._temperature, 

diff --git a/calphy/liquid.py b/calphy/liquid.py
@@ -70,7 +70,7 @@ def melt_structure(self, lmp):
 
             self.logger.info("Starting melting cycle with thigh temp %f, factor %f"%(self.calc._temperature_high, thmult))
             factor = (self.calc._temperature_high/self.calc._temperature)*thmult
-            lmp.velocity("all create", self.calc._temperature*factor, np.random.randint(0, 10000))
+            lmp.velocity("all create", self.calc._temperature*factor, np.random.randint(1, 10000))
             self.fix_nose_hoover(lmp, temp_start_factor=factor, temp_end_factor=factor)
             lmp.run(int(self.calc.md.n_small_steps))
             self.unfix_nose_hoover(lmp)
@@ -122,7 +122,7 @@ def run_averaging(self):
         lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
 
         #Melt regime for the liquid
-        lmp.velocity("all create", self.calc._temperature_high, np.random.randint(0, 10000))
+        lmp.velocity("all create", self.calc._temperature_high, np.random.randint(1, 10000))
 
         #add some computes
         lmp.command("variable         mvol equal vol")
@@ -192,7 +192,7 @@ def run_integration(self, iteration=1):
 
         lmp.command("fix              f1 all nve")
         lmp.command("fix              f2 all langevin %f %f %f %d zero yes"%(self.calc._temperature, self.calc._temperature, self.calc.md.thermostat_damping[1], 
-                                        np.random.randint(0, 10000)))
+                                        np.random.randint(1, 10000)))
         lmp.command("run               %d"%self.calc.n_equilibration_steps)
 
         lmp.command("unfix            f1")
@@ -225,11 +225,11 @@ def run_integration(self, iteration=1):
         lmp.command("thermo           1000")
 
 
-        lmp.command("velocity         all create %f %d mom yes rot yes dist gaussian"%(self.calc._temperature, np.random.randint(0, 10000)))
+        lmp.command("velocity         all create %f %d mom yes rot yes dist gaussian"%(self.calc._temperature, np.random.randint(1, 10000)))
 
         lmp.command("fix              f1 all nve")
         lmp.command("fix              f2 all langevin %f %f %f %d zero yes"%(self.calc._temperature, self.calc._temperature, self.calc.md.thermostat_damping[1], 
-                                        np.random.randint(0, 10000)))
+                                        np.random.randint(1, 10000)))
         lmp.command("compute          Tcm all temp/com")
         lmp.command("fix_modify       f2 temp Tcm")
 
@@ -254,7 +254,7 @@ def run_integration(self, iteration=1):
 
         lmp.command("fix              f1 all nve")
         lmp.command("fix              f2 all langevin %f %f %f %d zero yes"%(self.calc._temperature, self.calc._temperature, self.calc.md.thermostat_damping[1], 
-                                        np.random.randint(0, 10000)))
+                                        np.random.randint(1, 10000)))
         lmp.command("fix_modify       f2 temp Tcm")
 
         lmp.command("run               %d"%self.calc.n_equilibration_steps)
@@ -286,7 +286,7 @@ def run_integration(self, iteration=1):
 
         lmp.command("fix              f1 all nve")
         lmp.command("fix              f2 all langevin %f %f %f %d zero yes"%(self.calc._temperature, self.calc._temperature, self.calc.md.thermostat_damping[1], 
-                                        np.random.randint(0, 10000)))
+                                        np.random.randint(1, 10000)))
         lmp.command("fix_modify       f2 temp Tcm")
 
         lmp.command("fix              f3 all print 1 \"${dU1} ${dU2} ${flambda}\" screen no file backward_%d.dat"%iteration)

diff --git a/calphy/phase.py b/calphy/phase.py
@@ -57,7 +57,8 @@ def __init__(self, calculation=None, simfolder=None):
         self.logger = ph.prepare_log(logfile)
 
         self.logger.info("---------------input file----------------")
-        self.logger.info(yaml.safe_dump(self.calc.to_dict()))
+        self.logger.info("commented out as causes crash when we're expanding the T range after a fail run")
+       # self.logger.info(yaml.safe_dump(self.calc.to_dict()))
         self.logger.info("------------end of input file------------")
 
         if self.calc._pressure is None:
@@ -355,7 +356,7 @@ def run_zero_pressure_equilibration(self, lmp):
         Each method should close all the fixes. Run a small eqbr routine to achieve zero pressure        
         """
         #set velocity
-        lmp.command("velocity         all create %f %d"%(self.calc._temperature, np.random.randint(0, 10000)))
+        lmp.command("velocity         all create %f %d"%(self.calc._temperature, np.random.randint(1, 10000)))
 
         #apply fixes depending on thermostat/barostat
         if self.calc.equilibration_control == "nose-hoover":
@@ -397,7 +398,7 @@ def run_finite_pressure_equilibration(self, lmp):
         can prevent the issue.         
         """
         #create velocity
-        lmp.command("velocity         all create %f %d"%(0.25*self.calc._temperature, np.random.randint(0, 10000)))
+        lmp.command("velocity         all create %f %d"%(0.25*self.calc._temperature, np.random.randint(1, 10000)))
 
         #for Nose-Hoover thermo/baro combination
         if self.calc.equilibration_control == "nose-hoover":
@@ -578,7 +579,7 @@ def run_constrained_pressure_convergence(self, lmp):
     def run_iterative_constrained_pressure_convergence(self, lmp):
         """
         """
-        lmp.command("velocity         all create %f %d"%(self.calc._temperature, np.random.randint(0, 10000)))
+        lmp.command("velocity         all create %f %d"%(self.calc._temperature, np.random.randint(1, 10000)))
         lmp.command("fix              1 all nvt temp %f %f %f"%(self.calc._temperature, self.calc._temperature, self.calc.md.thermostat_damping[1]))
         lmp.command("thermo_style     custom step pe press vol etotal temp lx ly lz")
         lmp.command("thermo           10")
@@ -652,7 +653,7 @@ def finalise_pressure(self,):
     def run_minimal_constrained_pressure_convergence(self, lmp):
         """
         """
-        lmp.command("velocity         all create %f %d"%(self.calc._temperature, np.random.randint(0, 10000)))
+        lmp.command("velocity         all create %f %d"%(self.calc._temperature, np.random.randint(1, 10000)))
         lmp.command("fix              1 all nvt temp %f %f %f"%(self.calc._temperature, self.calc._temperature, self.calc.md.thermostat_damping[1]))
         lmp.command("thermo_style     custom step pe press vol etotal temp lx ly lz")
         lmp.command("thermo           10")
@@ -832,7 +833,7 @@ def reversible_scaling(self, iteration=1):
         lmp.command("thermo            10000")
 
         #create velocity and equilibriate
-        lmp.command("velocity          all create %f %d mom yes rot yes dist gaussian"%(t0, np.random.randint(0, 10000)))   
+        lmp.command("velocity          all create %f %d mom yes rot yes dist gaussian"%(t0, np.random.randint(1, 10000)))
         lmp.command("run               %d"%self.calc.n_equilibration_steps)
 
         lmp.command("variable         flambda equal ramp(${li},${lf})")
@@ -1131,4 +1132,4 @@ def integrate_pressure_scaling(self, return_values=False):
             nsims=self.calc.n_iterations, return_values=return_values)
 
         if return_values:
-            return res    
+    return res
diff --git a/calphy/solid.py b/calphy/solid.py
@@ -398,7 +398,7 @@ def run_integration(self, iteration=1):
 
         #apply temp fix
         lmp.command("fix               f3 all langevin %f %f %f %d zero yes"%(self.calc._temperature, self.calc._temperature, self.calc.md.thermostat_damping[1], 
-                                        np.random.randint(0, 10000)))
+                                        np.random.randint(1, 10000)))
 
         #compute com and apply to fix
         lmp.command("compute           Tcm all temp/com")
@@ -416,7 +416,7 @@ def run_integration(self, iteration=1):
         lmp.command("thermo            10000")
 
         #Create velocity
-        lmp.command("velocity          all create %f %d mom yes rot yes dist gaussian"%(self.calc._temperature, np.random.randint(0, 10000)))
+        lmp.command("velocity          all create %f %d mom yes rot yes dist gaussian"%(self.calc._temperature, np.random.randint(1, 10000)))
 
         #reapply 
         for i in range(self.calc.n_elements):