diff --git a/python/BioSimSpace/Align/_align.py b/python/BioSimSpace/Align/_align.py index bdbe025cb..3b9c6f28d 100644 --- a/python/BioSimSpace/Align/_align.py +++ b/python/BioSimSpace/Align/_align.py @@ -725,6 +725,9 @@ def matchAtoms( complete_rings_only=True, max_scoring_matches=1000, roi=None, + prune_perturbed_constraints=False, + prune_crossing_constraints=False, + prune_atom_types=False, property_map0={}, property_map1={}, ): @@ -787,6 +790,20 @@ def matchAtoms( The region of interest to match. Consists of a list of ROI residue indices. + prune_perturbed_constraints : bool + Whether to remove hydrogen atoms that are perturbed to heavy atoms + from the mapping. This option should be True when creating mappings + to use with the AMBER engine. + + prune_crossing_constraints : bool + Whether to remove atoms from the mapping such that there are no + constraints between dummy and non-dummy atoms. This option should + be True when creating mappings to use with the AMBER engine. + + prune_atom_types : bool + Whether to remove atoms from the mapping such that there are no + changing atom types. + property_map0 : dict A dictionary that maps "properties" in molecule0 to their user defined values. This allows the user to refer to properties @@ -859,6 +876,9 @@ def matchAtoms( timeout=timeout, complete_rings_only=complete_rings_only, max_scoring_matches=max_scoring_matches, + prune_perturbed_constraints=prune_perturbed_constraints, + prune_crossing_constraints=prune_crossing_constraints, + prune_atom_types=prune_atom_types, property_map0=property_map0, property_map1=property_map1, ) @@ -867,6 +887,9 @@ def matchAtoms( molecule0=molecule0, molecule1=molecule1, roi=roi, + prune_perturbed_constraints=prune_perturbed_constraints, + prune_crossing_constraints=prune_crossing_constraints, + prune_atom_types=prune_atom_types, use_kartograf=False, kartograf_kwargs={}, ) @@ -875,15 +898,19 @@ def matchAtoms( def _matchAtoms( molecule0, molecule1, - scoring_function, - matches, - return_scores, - prematch, - timeout, - complete_rings_only, - max_scoring_matches, - property_map0, - property_map1, + scoring_function="rmsd_align", + matches=1, + return_scores=False, + prematch={}, + timeout=5 * _Units.Time.second, + complete_rings_only=True, + max_scoring_matches=1000, + roi=None, + prune_perturbed_constraints=False, + prune_crossing_constraints=False, + prune_atom_types=False, + property_map0={}, + property_map1={}, ): # A list of supported scoring functions. scoring_functions = ["RMSD", "RMSDALIGN", "RMSDFLEXALIGN"] @@ -935,12 +962,24 @@ def _matchAtoms( if not isinstance(timeout, _Units.Time._Time): raise TypeError("'timeout' must be of type 'BioSimSpace.Types.Time'") + if not isinstance(complete_rings_only, bool): + raise TypeError("'complete_rings_only' must be of type 'bool'") + if not type(max_scoring_matches) is int: raise TypeError("'max_scoring_matches' must be of type 'int'") if max_scoring_matches <= 0: raise ValueError("'max_scoring_matches' must be >= 1.") + if not isinstance(prune_perturbed_constraints, bool): + raise TypeError("'prune_perturbed_constraints' must be of type 'bool'") + + if not isinstance(prune_crossing_constraints, bool): + raise TypeError("'prune_crossing_constraints' must be of type 'bool'") + + if not isinstance(prune_atom_types, bool): + raise TypeError("'prune_atom_types' must be of type 'bool'") + if not isinstance(property_map0, dict): raise TypeError("'property_map0' must be of type 'dict'") @@ -1069,6 +1108,22 @@ def _matchAtoms( property_map1, ) + # Optionally post-process the MCS for use with AMBER. + if prune_perturbed_constraints: + mappings = [ + _prune_perturbed_constraints(molecule0, molecule1, mapping) + for mapping in mappings + ] + if prune_crossing_constraints: + mappings = [ + _prune_crossing_constraints(molecule0, molecule1, mapping) + for mapping in mappings + ] + if prune_atom_types: + mappings = [ + _prune_atom_types(molecule0, molecule1, mapping) for mapping in mappings + ] + if matches == 1: if return_scores: return (mappings[0], scores[0]) @@ -1107,7 +1162,7 @@ def _kartograf_map(molecule0, molecule1, kartograf_kwargs): The kartograf mapping object. """ - # try to import kartograf + # Try to import kartograf. try: from kartograf.atom_aligner import align_mol_shape as _align_mol_shape from kartograf.atom_mapping_scorer import ( @@ -1138,22 +1193,21 @@ def _kartograf_map(molecule0, molecule1, kartograf_kwargs): rdkit_mol1 = _toRDKit(molecule1) rdkit_mols = [rdkit_mol0, rdkit_mol1] - # build small molecule components + # Build small molecule components. mol0, mol1 = [_SmallMoleculeComponent.from_rdkit(m) for m in rdkit_mols] - # align molecules first + # Align molecules first. a_mol1 = _align_mol_shape(mol1, ref_mol=mol0) - # build Kartograf Atom Mapper + # Build Kartograf Atom Mapper. mapper = KartografAtomMapper(**kartograf_kwargs) - # get the mapping + # Get the mapping. kartograf_mapping = next(mapper.suggest_mappings(mol0, a_mol1)) - # score the mapping + # Score the mapping. rmsd_scorer = _MappingRMSDScorer() score = rmsd_scorer(mapping=kartograf_mapping) - print(f"RMSD score: {score:.2f}") return kartograf_mapping @@ -1162,6 +1216,9 @@ def _roiMatch( molecule0, molecule1, roi, + prune_perturbed_constraints=False, + prune_crossing_constraints=False, + prune_atom_types=False, **kwargs, ): """ @@ -1184,6 +1241,20 @@ def _roiMatch( The region of interest to match. Consists of a list of ROI residue indices + prune_perturbed_constraints : bool + Whether to remove hydrogen atoms that are perturbed to heavy atoms + from the mapping. This option should be True when creating mappings + to use with the AMBER engine. + + prune_crossing_constraints : bool + Whether to remove atoms from the mapping such that there are no + constraints between dummy and non-dummy atoms. This option should + be True when creating mappings to use with the AMBER engine. + + prune_atom_types : bool + Whether to remove atoms from the mapping such that there are no + changing atom types. + use_kartograf : bool, optional, default=False If set to True, will use the kartograf algorithm to match the molecules. @@ -1255,7 +1326,16 @@ def _roiMatch( else: _validate_roi([molecule0, molecule1], roi) - # Check kwargs + if not isinstance(prune_perturbed_constraints, bool): + raise TypeError("'prune_perturbed_constraints' must be of type 'bool'") + + if not isinstance(prune_crossing_constraints, bool): + raise TypeError("'prune_crossing_constraints' must be of type 'bool'") + + if not isinstance(prune_atom_types, bool): + raise TypeError("'prune_atom_types' must be of type 'bool'") + + # Check kwargs. use_kartograf = kwargs.get("use_kartograf", False) kartograf_kwargs = kwargs.get("kartograf_kwargs", {}) @@ -1324,7 +1404,7 @@ def _roiMatch( res0_idx = [a.index() for a in molecule0_roi] res1_idx = [a.index() for a in molecule1_roi] - # Extract the residues of interest from the molecules + # Extract the residues of interest from the molecules. res0_extracted = molecule0.extract(res0_idx) res1_extracted = molecule1.extract(res1_idx) @@ -1414,13 +1494,24 @@ def _roiMatch( ) ) - # Combine the dictionaries to get the full mapping + # Combine the dictionaries to get the full mapping. full_mapping = { **pre_roi_mapping, **absolute_roi_mapping, **after_roi_mapping, } + # Optionally post-process the MCS for use with AMBER. + + if prune_perturbed_constraints: + full_mapping = _prune_perturbed_constraints(molecule0, molecule1, full_mapping) + + if prune_crossing_constraints: + full_mapping = _prune_crossing_constraints(molecule0, molecule1, full_mapping) + + if prune_atom_types: + full_mapping = _prune_atom_types(molecule0, molecule1, full_mapping) + return full_mapping @@ -2905,12 +2996,11 @@ def _validate_roi(molecules, roi): Parameters ---------- - molecules : list - Consits of a list of :class:`Molecule `. + molecules : list[:class:`Molecule `] + A list of molecules. roi : list - The region of interest. - Consists of a list of ROI residue indices. + A list of residue indices. """ if not isinstance(roi, list): @@ -2984,3 +3074,166 @@ def _from_sire_mapping(sire_mapping): mapping[idx0.value()] = idx1.value() return mapping + + +def _prune_perturbed_constraints(molecule0, molecule1, mapping): + """ + Prunes the maximum common substructure (MCS) so that no hydrogen + bond constraints are perturbed. + + Parameters + ---------- + + molecule0 : class:`Molecule + The first molecule (used at lambda = 0). + + molecule1 : class:`Molecule + The second molecule (used at lambda = 1). + + mapping : dict(int, int) + A maximum common substructure mapping between both molecules, as + generated by e.g. BioSimSpace.Align.matchAtoms(). + + Returns + ------- + + new_mapping : dict(int, int) + The pruned MCS. + """ + new_mapping = {} + + # Store a hydrogen element. + hydrogen = _SireMol.Element("H") + + for idx0, idx1 in mapping.items(): + atom0 = molecule0.getAtoms()[idx0] + atom1 = molecule1.getAtoms()[idx1] + elem0 = atom0._sire_object.property("element") + elem1 = atom1._sire_object.property("element") + elems = {elem0, elem1} + + # Make sure we are not matching a hydrogen to a non-hydrogen. + if not (hydrogen in elems and len(elems) > 1): + new_mapping[idx0] = idx1 + + return new_mapping + + +def _prune_crossing_constraints(molecule0, molecule1, mapping): + """ + Prunes the maximum common substructure (MCS) mapping so that there are no + constrained bonds between a common core and a softcore atom. + + Parameters + ---------- + + molecule0 : class:`Molecule + The first molecule (used at lambda = 0). + + molecule1 : class:`Molecule + The second molecule (used at lambda = 1). + + mapping : dict(int, int) + A maximum common substructure mapping between both molecules, as + generated by e.g. BioSimSpace.Align.matchAtoms(). + + Returns + ------- + + new_mapping : dict(int, int) + The pruned mapping. + """ + + # Get the connectivity of the molecules. + connectivity0 = _SireMol.Connectivity( + molecule0._sire_object, _SireMol.CovalentBondHunter() + ) + connectivity1 = _SireMol.Connectivity( + molecule1._sire_object, _SireMol.CovalentBondHunter() + ) + + # Store a hydrogen element. + hydrogen = _SireMol.Element("H") + + while True: + new_mapping = {} + + for idx0, idx1 in mapping.items(): + # Get the relevant atom and whether it's a hydrogen. + atom0 = molecule0._sire_object.atom(_SireMol.AtomIdx(idx0)) + atom1 = molecule1._sire_object.atom(_SireMol.AtomIdx(idx1)) + is_H0 = atom0.property("element") == hydrogen + is_H1 = atom1.property("element") == hydrogen + + # Get the neighbours to the atom + neighbours0 = [ + molecule0._sire_object.atom(i) + for i in connectivity0.connectionsTo(_SireMol.AtomIdx(idx0)) + ] + neighbours1 = [ + molecule1._sire_object.atom(i) + for i in connectivity1.connectionsTo(_SireMol.AtomIdx(idx1)) + ] + + # Determine whether there are any constrained bonds between the + # MCS and softcore part. + any_Hdummies0 = any( + (atom.property("element") == hydrogen or is_H0) + and atom.index().value() not in mapping.keys() + for atom in neighbours0 + ) + any_Hdummies1 = any( + (atom.property("element") == hydrogen or is_H1) + and atom.index().value() not in mapping.values() + for atom in neighbours1 + ) + + if not any_Hdummies0 and not any_Hdummies1: + new_mapping[idx0] = idx1 + + # We stop iterating if the pruned mapping is the same as the input one. + if new_mapping == mapping: + mapping = new_mapping + break + else: + mapping = new_mapping + + return mapping + + +def _prune_atom_types(molecule0, molecule1, mapping): + """ + Prunes the maximum common substructure (MCS) so that there are no + atoms changing type. + + Parameters + ---------- + + molecule0 : class:`Molecule + The first molecule (used at lambda = 0). + + molecule1 : class:`Molecule + The second molecule (used at lambda = 1). + + mapping : dict(int, int) + A maximum common substructure mapping between both molecules, as + generated by e.g. BioSimSpace.Align.matchAtoms(). + + Returns + ------- + + new_mapping : dict(int, int) + The pruned mapping. + """ + new_mapping = {} + + for idx0, idx1 in mapping.items(): + atom0 = molecule0.getAtoms()[idx0] + atom1 = molecule1.getAtoms()[idx1] + elem0 = atom0._sire_object.property("element") + elem1 = atom1._sire_object.property("element") + + if elem0 == elem1: + new_mapping[idx0] = idx1 + + return new_mapping diff --git a/python/BioSimSpace/FreeEnergy/_relative.py b/python/BioSimSpace/FreeEnergy/_relative.py index 361b2b40a..e41027ca2 100644 --- a/python/BioSimSpace/FreeEnergy/_relative.py +++ b/python/BioSimSpace/FreeEnergy/_relative.py @@ -163,7 +163,7 @@ def __init__( engine : str The molecular dynamics engine used to run the simulation. Available - options are "GROMACS", or "SOMD". If this argument is omitted then + options are "AMBER", "GROMACS", or "SOMD". If this argument is omitted then BioSimSpace will choose an appropriate engine for you. setup_only : bool @@ -1158,8 +1158,8 @@ def _preprocess_data(data, estimator, **kwargs): # Assign defaults in case not passed via kwargs. auto_eq = False stat_ineff = False - truncate = False - truncate_keep = "start" + truncate_upper = 100 + truncate_lower = 0 # Parse kwargs. for key, value in kwargs.items(): @@ -1168,35 +1168,31 @@ def _preprocess_data(data, estimator, **kwargs): auto_eq = value if key == "STATISTICALINEFFICIENCY": stat_ineff = value - if key == "TRUNCATEPERCENTAGE": - truncate = value - if key == "TRUNCATEKEEP": - truncate_keep = value + if key == "TRUNCATEUPPER": + truncate_upper = value + if key == "TRUNCATELOWER": + truncate_lower = value - # First truncate data. + # Copy the data. raw_data = data - if truncate: - # Get the upper and lower bounds for truncate. - data_len = len(data[0]) - data_step = round((data[0].index[-1][0] - data[0].index[-2][0]), 1) - data_kept = data_len * (truncate / 100) - data_time = data_kept * data_step - if truncate_keep == "start": - truncate_lower = 0 - truncate_upper = data_time - data_step - if truncate_keep == "end": - truncate_lower = (data_len * data_step) - data_time - truncate_upper = (data_len * data_step) - data_step - try: - data = [ - _slicing(i, lower=truncate_lower, upper=truncate_upper) - for i in raw_data - ] - except: - _warnings.warn("Could not truncate data.") - data = raw_data - else: + # Data length. + data_len = len(data[0]) + + # Step size. + data_step = round((data[0].index[-1][0] - data[0].index[-2][0]), 1) + + # Get the upper and lower bounds for truncate. + truncate_lower = (data_len * (truncate_lower / 100)) * data_step + truncate_upper = (data_len * (truncate_upper / 100)) * data_step + + try: + data = [ + _slicing(i, lower=truncate_lower, upper=truncate_upper) + for i in raw_data + ] + except: + _warnings.warn("Could not truncate data.") data = raw_data # The decorrelate function calls either autoequilibration or statistical_inefficiency @@ -1230,10 +1226,7 @@ def _preprocess_data(data, estimator, **kwargs): ) sampled_data = data - # Concatanate in alchemlyb format. - processed_data = _alchemlyb.concat(sampled_data) - - return processed_data + return sampled_data @staticmethod def _analyse_internal(files, temperatures, lambdas, engine, estimator, **kwargs): @@ -1319,6 +1312,7 @@ def _analyse_internal(files, temperatures, lambdas, engine, estimator, **kwargs) # Preprocess the data. try: processed_data = Relative._preprocess_data(data, estimator, **kwargs) + processed_data = _alchemlyb.concat(processed_data) except: _warnings.warn("Could not preprocess the data!") processed_data = _alchemlyb.concat(data) diff --git a/python/BioSimSpace/Process/_amber.py b/python/BioSimSpace/Process/_amber.py index 5c3a5b602..1bf0aebd9 100644 --- a/python/BioSimSpace/Process/_amber.py +++ b/python/BioSimSpace/Process/_amber.py @@ -277,6 +277,14 @@ def _setup(self, **kwargs): "perturbable molecule!" ) + # Make sure the protocol is valid. + if self._protocol.getPerturbationType() != "full": + raise NotImplementedError( + "AMBER currently only supports the 'full' perturbation " + "type. Please use engine='SOMD' when running multistep " + "perturbation types." + ) + # If this is vacuum simulation with pmemd.cuda then # we need to add a simulation box. if self._is_vacuum and self._is_pmemd_cuda: diff --git a/python/BioSimSpace/Process/_openmm.py b/python/BioSimSpace/Process/_openmm.py index 327c8e2cb..0a984af33 100644 --- a/python/BioSimSpace/Process/_openmm.py +++ b/python/BioSimSpace/Process/_openmm.py @@ -1488,13 +1488,12 @@ def getFrame(self, index): if not type(index) is int: raise TypeError("'index' must be of type 'int'") - max_index = ( - int( - (self._protocol.getRunTime() / self._protocol.getTimeStep()) - / self._protocol.getRestartInterval() - ) - - 1 + + max_index = int( + (self._protocol.getRunTime() / self._protocol.getTimeStep()) + / self._protocol.getRestartInterval() ) + if index < 0 or index > max_index: raise ValueError(f"'index' must be in range [0, {max_index}].") diff --git a/python/BioSimSpace/Protocol/_equilibration.py b/python/BioSimSpace/Protocol/_equilibration.py index 5bed6f73b..aac74933a 100644 --- a/python/BioSimSpace/Protocol/_equilibration.py +++ b/python/BioSimSpace/Protocol/_equilibration.py @@ -42,14 +42,15 @@ class Equilibration(_Protocol, _PositionRestraintMixin): def __init__( self, timestep=_Types.Time(2, "femtosecond"), - runtime=_Types.Time(0.2, "nanoseconds"), + runtime=_Types.Time(1, "nanoseconds"), temperature_start=_Types.Temperature(300, "kelvin"), temperature_end=_Types.Temperature(300, "kelvin"), temperature=None, pressure=None, thermostat_time_constant=_Types.Time(1, "picosecond"), - report_interval=100, - restart_interval=500, + report_interval=200, + restart_interval=1000, + restart=False, restraint=None, force_constant=10 * _Units.Energy.kcal_per_mol / _Units.Area.angstrom2, ): @@ -89,6 +90,9 @@ def __init__( The frequency at which restart configurations and trajectory frames are saved. (In integration steps.) + restart : bool + Whether this is a continuation of a previous simulation. + restraint : str, [int] The type of restraint to perform. This should be one of the following options: @@ -155,6 +159,9 @@ def __init__( # Set the restart interval. self.setRestartInterval(restart_interval) + # Set the restart flag. + self.setRestart(restart) + # Set the posistion restraint. _PositionRestraintMixin.__init__(self, restraint, force_constant) @@ -169,7 +176,7 @@ def _get_parm(self): f"thermostat_time_constant={self._thermostat_time_constant}, " f"report_interval={self._report_interval}, " f"restart_interval={self._restart_interval}, " - + _PositionRestraintMixin._get_parm(self) + f"restart={self._restart}, " + _PositionRestraintMixin._get_parm(self) ) def __str__(self): @@ -204,6 +211,7 @@ def __eq__(self, other): and self._thermostat_time_constant == other._thermostat_time_constant and self._report_interval == other._report_interval and self._restart_interval == other._restart_interval + and self._restart == other._restart and _PositionRestraintMixin.__eq__(self, other) ) @@ -484,6 +492,34 @@ def setRestartInterval(self, restart_interval): self._restart_interval = restart_interval + def isRestart(self): + """ + Return whether this restart simulation. + + Returns + ------- + + is_restart : bool + Whether this is a restart simulation. + """ + return self._restart + + def setRestart(self, restart): + """ + Set the restart flag. + + Parameters + ---------- + + restart : bool + Whether this is a restart simulation. + """ + if isinstance(restart, bool): + self._restart = restart + else: + _warnings.warn("Non-boolean restart flag. Defaulting to False!") + self._restart = False + def isConstantTemp(self): """ Return whether the protocol has a constant temperature. diff --git a/python/BioSimSpace/Protocol/_free_energy_equilibration.py b/python/BioSimSpace/Protocol/_free_energy_equilibration.py index 1ab69e2e2..57fdad8dc 100644 --- a/python/BioSimSpace/Protocol/_free_energy_equilibration.py +++ b/python/BioSimSpace/Protocol/_free_energy_equilibration.py @@ -44,7 +44,7 @@ def __init__( max_lam=1.0, num_lam=11, timestep=_Types.Time(2, "femtosecond"), - runtime=_Types.Time(0.2, "nanoseconds"), + runtime=_Types.Time(1, "nanoseconds"), temperature_start=_Types.Temperature(300, "kelvin"), temperature_end=_Types.Temperature(300, "kelvin"), temperature=None, @@ -52,6 +52,7 @@ def __init__( thermostat_time_constant=_Types.Time(1, "picosecond"), report_interval=200, restart_interval=1000, + restart=False, restraint=None, force_constant=10 * _Units.Energy.kcal_per_mol / _Units.Area.angstrom2, perturbation_type="full", @@ -106,6 +107,9 @@ def __init__( The frequency at which restart configurations and trajectory frames are saved. (In integration steps.) + restart : bool + Whether this is a continuation of a previous simulation. + restraint : str, [int] The type of restraint to perform. This should be one of the following options: @@ -217,4 +221,5 @@ def _to_regular_protocol(self): restart_interval=self.getRestartInterval(), restraint=self.getRestraint(), force_constant=self.getForceConstant(), + restart=self.isRestart(), ) diff --git a/python/BioSimSpace/Protocol/_production.py b/python/BioSimSpace/Protocol/_production.py index 481bfc1e4..811522cfd 100644 --- a/python/BioSimSpace/Protocol/_production.py +++ b/python/BioSimSpace/Protocol/_production.py @@ -46,8 +46,8 @@ def __init__( temperature=_Types.Temperature(300, "kelvin"), pressure=_Types.Pressure(1, "atmosphere"), thermostat_time_constant=_Types.Time(1, "picosecond"), - report_interval=100, - restart_interval=100, + report_interval=200, + restart_interval=1000, restart=False, restraint=None, force_constant=10 * _Units.Energy.kcal_per_mol / _Units.Area.angstrom2, diff --git a/python/BioSimSpace/_Config/_config.py b/python/BioSimSpace/_Config/_config.py index 9005dbb78..420c9a9cd 100644 --- a/python/BioSimSpace/_Config/_config.py +++ b/python/BioSimSpace/_Config/_config.py @@ -213,6 +213,6 @@ def steps(self): steps = self._protocol.getSteps() else: steps = _math.ceil( - self._protocol.getRunTime() / self._protocol.getTimeStep() + round(self._protocol.getRunTime() / self._protocol.getTimeStep(), 6) ) return steps diff --git a/python/BioSimSpace/_Config/_gromacs.py b/python/BioSimSpace/_Config/_gromacs.py index 1355b274c..dae45c24b 100644 --- a/python/BioSimSpace/_Config/_gromacs.py +++ b/python/BioSimSpace/_Config/_gromacs.py @@ -100,23 +100,12 @@ def createConfig(self, version=None, extra_options={}, extra_lines=[]): if not all(isinstance(line, str) for line in extra_lines): raise TypeError("Lines in 'extra_lines' must be of type 'str'.") - # Make sure the report interval is a multiple of nstcalcenergy. - if isinstance(self._protocol, _FreeEnergyMixin): - nstcalcenergy = 250 - else: - nstcalcenergy = 100 - report_interval = self.reportInterval() - if report_interval % nstcalcenergy != 0: - report_interval = nstcalcenergy * _math.ceil( - report_interval / nstcalcenergy - ) - # Define some miscellaneous defaults. protocol_dict = { # Interval between writing to the log file. - "nstlog": report_interval, + "nstlog": self.reportInterval(), # Interval between writing to the energy file. - "nstenergy": report_interval, + "nstenergy": self.reportInterval(), # Interval between writing to the trajectory file. "nstxout-compressed": self.restartInterval(), } @@ -124,6 +113,8 @@ def createConfig(self, version=None, extra_options={}, extra_lines=[]): # Minimisation. if isinstance(self._protocol, _Protocol.Minimisation): protocol_dict["integrator"] = "steep" + # Maximum step size in nanometers. + protocol_dict["emstep"] = "0.001" else: # Timestep in picoseconds timestep = self._protocol.getTimeStep().picoseconds().value() @@ -212,7 +203,7 @@ def createConfig(self, version=None, extra_options={}, extra_lines=[]): # Temperature control. if not isinstance(self._protocol, _Protocol.Minimisation): if isinstance(self._protocol, _FreeEnergyMixin): - # Langevin dynamics. + # Leap-frog stochastic dynamics integrator. protocol_dict["integrator"] = "sd" else: # Leap-frog molecular dynamics. @@ -260,6 +251,24 @@ def createConfig(self, version=None, extra_options={}, extra_lines=[]): "%.2f" % self._protocol.getTemperature().kelvin().value() ) + # Set as a continuation run. + if self.isRestart(): + protocol_dict["continuation"] = "yes" + protocol_dict["gen-vel"] = "no" + # Generate velocities. + else: + protocol_dict["continuation"] = "no" + protocol_dict["gen-vel"] = "yes" + protocol_dict["gen-seed"] = "-1" + if isinstance(self._protocol, _Protocol.Equilibration): + protocol_dict["gen-temp"] = ( + "%.2f" % self._protocol.getStartTemperature().kelvin().value() + ) + else: + protocol_dict["gen-temp"] = ( + "%.2f" % self._protocol.getTemperature().kelvin().value() + ) + # Free energies. if isinstance(self._protocol, _FreeEnergyMixin): # Extract the lambda array. @@ -280,10 +289,14 @@ def createConfig(self, version=None, extra_options={}, extra_lines=[]): protocol_dict["couple-lambda1"] = "vdw-q" # Write all lambda values. protocol_dict["calc-lambda-neighbors"] = -1 - # Calculate energies every 250 steps. - protocol_dict["nstcalcenergy"] = 250 - # Write gradients every 250 steps. - protocol_dict["nstdhdl"] = 250 + # Calculate energies at the report interval. + protocol_dict["nstcalcenergy"] = self.reportInterval() + # Write gradients at the report interval. + protocol_dict["nstdhdl"] = self.reportInterval() + # Soft-core parameters. + protocol_dict["sc-alpha"] = "0.30" + protocol_dict["sc-sigma"] = "0.25" + protocol_dict["sc-coul"] = "yes" # Put everything together in a line-by-line format. total_dict = {**protocol_dict, **extra_options} diff --git a/python/BioSimSpace/_Config/_somd.py b/python/BioSimSpace/_Config/_somd.py index a5ba74106..6d4bc4118 100644 --- a/python/BioSimSpace/_Config/_somd.py +++ b/python/BioSimSpace/_Config/_somd.py @@ -30,6 +30,7 @@ import warnings as _warnings from .. import Protocol as _Protocol +from ..Protocol._free_energy_mixin import _FreeEnergyMixin from ..Protocol._position_restraint_mixin import _PositionRestraintMixin from ._config import Config as _Config @@ -112,51 +113,51 @@ def createConfig(self, extra_options={}, extra_lines=[]): report_interval = self._protocol.getReportInterval() restart_interval = self._protocol.getRestartInterval() - # Work out the number of cycles. - ncycles = ( - self._protocol.getRunTime() / self._protocol.getTimeStep() - ) / 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( - (self._protocol.getRunTime() / self._protocol.getTimeStep()) - / 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.FreeEnergyProduction): - 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. - if isinstance(self._protocol, _Protocol.FreeEnergyProduction): + # is report interval. Round down to the closest multiple. + if isinstance(self._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. @@ -250,8 +251,8 @@ def createConfig(self, extra_options={}, extra_lines=[]): protocol_dict["minimise"] = True # Handle hydrogen perturbations. protocol_dict["constraint"] = "hbonds-notperturbed" - # Write gradients every 250 steps. - protocol_dict["energy frequency"] = 250 + # Write gradients at the report interval. + protocol_dict["energy frequency"] = report_interval protocol = [str(x) for x in self._protocol.getLambdaValues()] protocol_dict["lambda array"] = ", ".join(protocol) diff --git a/tests/FreeEnergy/test_relative.py b/tests/FreeEnergy/test_relative.py index 6808b0ed7..aba68d96f 100644 --- a/tests/FreeEnergy/test_relative.py +++ b/tests/FreeEnergy/test_relative.py @@ -43,7 +43,7 @@ def expected_results(): """A dictionary of expected FEP results.""" return { - "amber": {"mbar": -12.5939, "ti": -13.6850}, + "amber": {"mbar": -12.7176, "ti": -13.8771}, "gromacs": {"mbar": -6.0238, "ti": -8.4158}, "somd": {"mbar": -6.3519, "ti": -6.3209}, } diff --git a/tests/Process/test_somd.py b/tests/Process/test_somd.py index f7eedc52c..0f077ca17 100644 --- a/tests/Process/test_somd.py +++ b/tests/Process/test_somd.py @@ -43,9 +43,7 @@ def test_free_energy(perturbable_system): """Test a free energy perturbation protocol.""" # Create a short FEP protocol. - protocol = BSS.Protocol.FreeEnergy( - runtime=0.1 * BSS.Units.Time.picosecond, report_interval=50, restart_interval=50 - ) + protocol = BSS.Protocol.FreeEnergy(runtime=0.1 * BSS.Units.Time.picosecond) # Run the process, check that it finished without error, and returns a system. run_process(perturbable_system, protocol)