From 8a323fcab26bc60802d921281b312c1c6a86c5fe Mon Sep 17 00:00:00 2001 From: finlayclark Date: Wed, 12 Feb 2025 21:22:22 +0000 Subject: [PATCH] Update Exscientia Sandpit SOMD config settings Update settings to match those inn the main version of the code: see https://github.com/OpenBioSim/biosimspace/commit/976ab368afcb8dc541f5eabc8fe7b7dbda6604a1 --- .../Sandpit/Exscientia/Protocol/_config.py | 55 ++++++++++--------- .../Exscientia/Protocol/test_config.py | 9 ++- 2 files changed, 36 insertions(+), 28 deletions(-) diff --git a/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py b/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py index bda35fb8d..0706510ba 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py +++ b/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()] diff --git a/tests/Sandpit/Exscientia/Protocol/test_config.py b/tests/Sandpit/Exscientia/Protocol/test_config.py index 32446dcda..1b97a5c63 100644 --- a/tests/Sandpit/Exscientia/Protocol/test_config.py +++ b/tests/Sandpit/Exscientia/Protocol/test_config.py @@ -480,7 +480,9 @@ 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 ) @@ -488,6 +490,11 @@ def test_turn_on_restraint_boresch(self, system_and_boresch_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"'