Fixed somd-freenrg OpenMM random number generator seeding bug.

CHANGELOG

@@ -3,7 +3,7 @@ SIRE changelog
 
 devel branch:
 
-    [2021.1.0] Feb 2021 - Added support for multiple combining rules in SOMD
+    [2021.1.0] Aug 2021 - Added support for multiple combining rules in SOMD
                (@SofiaBariami). Added support for triclinic simulation boxes.
                Convert Ryckaert-Bellememans form dihedral functions from
                GROMACS to Fourier series to allow conversion to AMBER format.
@@ -24,6 +24,9 @@ devel branch:
                conda version of libcpuid. Added support for generating PDB
                CONECT records from a Sire.Mol.Connectivity object. Fixed
                issue with PMEMD skipping torsions with zero periodicity.
+               Fixed random number seeding bug in somd-freenrg, which
+               resulted in the OpenMM generator being seeded with the same
+               seed for each cycle of the simulation.
 
     [2020.1.0] July 2020 - Fixed bug in WaterView program to ensure that a
                molecule is extracted from the returned list. Stable sorting

wrapper/Tools/OpenMMMD.py

@@ -1,4 +1,3 @@
-
 ####################################################################################################
 #                                                                                                  #
 #   RUN SCRIPT to perform an MD simulation in Sire with OpenMM                                     #
@@ -87,8 +86,10 @@
 
 nmoves = Parameter("nmoves", 1000, """Number of Molecular Dynamics moves to be performed during the simulation.""")
 
-random_seed = Parameter("random seed", None, """Random number seed. Set this if you
-                         want to have reproducible simulations.""")
+debug_seed = Parameter("debug seed", 0, """Debugging seed number seed. Set this if you
+                         want to reproduce a single cycle. Don't use this seed for production simulations
+                         since the same seed will be used for all cycles! A value of zero means that a unique
+                         seed will be generated for each cycle.""")
 
 ncycles = Parameter("ncycles", 1,
                     """The number of MD cycles. The total elapsed time will be nmoves*ncycles*timestep""")
@@ -455,7 +456,7 @@ def setupForcefields(system, space):
     return system
 
 
-def setupMoves(system, random_seed, GPUS):
+def setupMoves(system, debug_seed, GPUS):
 
     print("Setting up moves...")
 
@@ -507,11 +508,11 @@ def setupMoves(system, random_seed, GPUS):
     moves = WeightedMoves()
     moves.add(mdmove, 1)
 
-    if (not random_seed):
-        random_seed = RanGenerator().randInt(100000, 1000000)
-    print("Generated random seed number %d " % random_seed)
+    if debug_seed != 0:
+        debug_seed = RanGenerator().randInt(100000, 1000000)
+        print("Generated debugging seed number %d " % debug_seed)
 
-    moves.setGenerator(RanGenerator(random_seed))
+    moves.setGenerator(RanGenerator(debug_seed))
 
     return moves
 
@@ -1154,7 +1155,7 @@ def setupForceFieldsFreeEnergy(system, space):
     return system
 
 
-def setupMovesFreeEnergy(system, random_seed, GPUS, lam_val):
+def setupMovesFreeEnergy(system, debug_seed, GPUS, lam_val):
 
     print ("Setting up moves...")
 
@@ -1165,7 +1166,7 @@ def setupMovesFreeEnergy(system, random_seed, GPUS, lam_val):
     solute_fromdummy = system[MGName("solute_ref_fromdummy")]
 
     Integrator_OpenMM = OpenMMFrEnergyST(molecules, solute, solute_hard, solute_todummy, solute_fromdummy)
-    Integrator_OpenMM.setRandomSeed(random_seed)
+    Integrator_OpenMM.setRandomSeed(debug_seed)
     Integrator_OpenMM.setIntegrator(integrator_type.val)
     Integrator_OpenMM.setFriction(inverse_friction.val)  # Only meaningful for Langevin/Brownian integrators
     Integrator_OpenMM.setPlatform(platform.val)
@@ -1205,7 +1206,7 @@ def setupMovesFreeEnergy(system, random_seed, GPUS, lam_val):
     #This calls the OpenMMFrEnergyST initialise function
     Integrator_OpenMM.initialise()
     velocity_generator = MaxwellBoltzmann(temperature.val)
-    velocity_generator.setGenerator(RanGenerator(random_seed))
+    velocity_generator.setGenerator(RanGenerator(debug_seed))
 
     mdmove = MolecularDynamics(molecules, Integrator_OpenMM, timestep.val,
                               {"velocity generator":velocity_generator})
@@ -1215,12 +1216,11 @@ def setupMovesFreeEnergy(system, random_seed, GPUS, lam_val):
     moves = WeightedMoves()
     moves.add(mdmove, 1)
 
-    if (not random_seed):
-        random_seed = RanGenerator().randInt(100000, 1000000)
-
-    print("Generated random seed number %d " % random_seed)
+    if debug_seed != 0:
+        debug_seed = RanGenerator().randInt(100000, 1000000)
+        print("Generated debugging seed number %d " % debug_seed)
 
-    moves.setGenerator(RanGenerator(random_seed))
+    moves.setGenerator(RanGenerator(debug_seed))
 
     return moves
 
@@ -1429,12 +1429,11 @@ def run():
 
         system = setupForcefields(system, space)
 
-        if random_seed.val:
-            ranseed = random_seed.val
+        if debug_seed.val == 0:
+            ranseed = debug_seed.val
         else:
             ranseed = RanGenerator().randInt(100000, 1000000)
-
-        print("Setting up the simulation with random seed %s" % ranseed)
+            print("Setting up the simulation with debugging seed %s" % ranseed)
 
         moves = setupMoves(system, ranseed, gpu.val)
 
@@ -1587,12 +1586,11 @@ def runFreeNrg():
 
         system = setupForceFieldsFreeEnergy(system, space)
 
-        if random_seed.val:
-            ranseed = random_seed.val
+        if debug_seed.val == 0:
+            ranseed = debug_seed.val
         else:
             ranseed = RanGenerator().randInt(100000, 1000000)
-
-        print("Setting up the simulation with random seed %s" % ranseed)
+            print("Setting up the simulation with debugging seed %s" % ranseed)
 
         moves = setupMovesFreeEnergy(system, ranseed, gpu.val, lambda_val.val)