Merge pull request #390 from fjclark/devel

Update Exscientia Sandpit SOMD config settings

Author: lohedges

python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py

@@ -743,51 +743,52 @@ def generateSomdConfig(
             # Get the report and restart intervals.
             report_interval = self._report_interval
             restart_interval = self._restart_interval
-            runtime = self.protocol.getRunTime()
 
-            # Get the timestep.
-            timestep = self.protocol._timestep
-
-            # Work out the number of cycles.
-            ncycles = (runtime / timestep) / report_interval
-
-            # If the number of cycles isn't integer valued, adjust the report
-            # interval so that we match specified the run time.
-            if ncycles - _math.floor(ncycles) != 0:
-                ncycles = _math.floor(ncycles)
-                if ncycles == 0:
-                    ncycles = 1
-                report_interval = _math.ceil((runtime / timestep) / ncycles)
+            # Work out the number of cycles. We want one per nanosecond of simulation.
+            if self.protocol.getRunTime().nanoseconds().value() < 2:
+                ncycles = 1
+                moves_per_cycle = self._steps
+            else:
+                ncycles = _math.ceil(self.protocol.getRunTime().nanoseconds().value())
+                moves_per_cycle = _math.ceil(self._steps / ncycles)
 
-            # For free energy simulations, the report interval must be a multiple
-            # of the energy frequency which is 250 steps.
-            if isinstance(self.protocol, _Protocol._FreeEnergyMixin):
-                if report_interval % 250 != 0:
-                    report_interval = 250 * _math.ceil(report_interval / 250)
+            # The number of moves needs to be multiple of the report interval.
+            if moves_per_cycle % report_interval != 0:
+                moves_per_cycle = (
+                    _math.ceil(moves_per_cycle / report_interval) * report_interval
+                )
 
             # Work out the number of cycles per frame.
-            cycles_per_frame = restart_interval / report_interval
+            cycles_per_frame = restart_interval / moves_per_cycle
 
             # Work out whether we need to adjust the buffer frequency.
             buffer_freq = 0
             if cycles_per_frame < 1:
-                buffer_freq = cycles_per_frame * restart_interval
+                buffer_freq = restart_interval
                 cycles_per_frame = 1
-                self._buffer_freq = buffer_freq
             else:
                 cycles_per_frame = _math.floor(cycles_per_frame)
 
+            # Make sure that we aren't buffering more than 1000 frames per cycle.
+            if buffer_freq > 0 and moves_per_cycle / buffer_freq > 1000:
+                _warnings.warn(
+                    "Trajectory buffering will exceed limit. Reducing buffer frequency."
+                )
+                buffer_freq = moves_per_cycle / 1000
+
             # For free energy simulations, the buffer frequency must be an integer
             # multiple of the frequency at which free energies are written, which
-            # is 250 steps. Round down to the closest multiple.
+            # is report interval. Round down to the closest multiple.
             if isinstance(self.protocol, _Protocol._FreeEnergyMixin):
                 if buffer_freq > 0:
-                    buffer_freq = 250 * _math.floor(buffer_freq / 250)
+                    buffer_freq = report_interval * _math.floor(
+                        buffer_freq / report_interval
+                    )
 
             # The number of SOMD cycles.
-            protocol_dict["ncycles"] = int(ncycles)
+            protocol_dict["ncycles"] = ncycles
             # The number of moves per cycle.
-            protocol_dict["nmoves"] = report_interval
+            protocol_dict["nmoves"] = moves_per_cycle
             # Cycles per trajectory write.
             protocol_dict["ncycles_per_snap"] = cycles_per_frame
             # Buffering frequency.
@@ -862,7 +863,7 @@ def generateSomdConfig(
                     "hbonds-notperturbed"  # Handle hydrogen perturbations.
                 )
                 protocol_dict["energy frequency"] = (
-                    250  # Write gradients every 250 steps.
+                    report_interval  # Write gradients at report interval.
                 )
 
             protocol = [str(x) for x in self.protocol.getLambdaValues()]

tests/Sandpit/Exscientia/Protocol/test_config.py

@@ -480,14 +480,21 @@ def system_and_boresch_restraint():
     def test_turn_on_restraint_boresch(self, system_and_boresch_restraint):
         """Test for turning on multiple distance restraints"""
         system, restraint = system_and_boresch_restraint
-        protocol = FreeEnergy(perturbation_type="restraint")
+        protocol = FreeEnergy(perturbation_type="restraint", 
+                              runtime=1*BSS.Units.Time.nanosecond,
+                              timestep=2*BSS.Units.Time.femtosecond)
         freenrg = BSS.FreeEnergy.AlchemicalFreeEnergy(
             system, protocol, engine="SOMD", restraint=restraint
         )
 
         # Test .cfg file
         with open(f"{freenrg._work_dir}/lambda_0.0000/somd.cfg", "r") as f:
             cfg_text = f.read()
+            assert "nmoves = 500000" in cfg_text
+            assert "ncycles_per_snap = 1" in cfg_text
+            assert "buffered coordinates frequency = 1000" in cfg_text
+            assert "timestep = 2.00 femtosecond" in cfg_text
+            assert "energy frequency = 200" in cfg_text
             assert "use boresch restraints = True" in cfg_text
             assert 'boresch restraints dictionary = {"anchor_points":{"r1":1,'
             ' "r2":2, "r3":3, "l1":4, "l2":5, "l3":6}, "equilibrium_values"'