Source code for BioSimSpace.Process._openmm

######################################################################
# BioSimSpace: Making biomolecular simulation a breeze!
#
# Copyright: 2017-2024
#
# Authors: Lester Hedges <lester.hedges@gmail.com>
#
# BioSimSpace is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with BioSimSpace. If not, see <http://www.gnu.org/licenses/>.
#####################################################################

"""Functionality for running simulations with OpenMM."""

__author__ = "Lester Hedges"
__email__ = "lester.hedges@gmail.com"

__all__ = ["OpenMM"]

from .._Utils import _try_import

import math as _math
import os as _os

_pygtail = _try_import("pygtail")
import sys as _sys
import shutil as _shutil
import timeit as _timeit
import warnings as _warnings

from sire.legacy import Base as _SireBase
from sire.legacy import IO as _SireIO
from sire.legacy import Units as _SireUnits


from .. import _isVerbose
from .._Exceptions import IncompatibleError as _IncompatibleError
from .._Exceptions import MissingSoftwareError as _MissingSoftwareError
from .._SireWrappers import System as _System
from ..Metadynamics import CollectiveVariable as _CollectiveVariable
from ..Protocol._position_restraint_mixin import _PositionRestraintMixin
from ..Types._type import Type as _Type
from .. import IO as _IO
from .. import Protocol as _Protocol
from .. import Trajectory as _Trajectory
from .. import Types as _Types
from .. import Units as _Units
from .. import _Utils

from . import _process

from ._plumed import Plumed as _Plumed


[docs] class OpenMM(_process.Process): """A class for running simulations using OpenMM.""" # Dictionary of platforms and their OpenMM keyword. _platforms = {"CPU": "CPU", "CUDA": "CUDA", "OPENCL": "OpenCL"} # Special cases for generate config when using ATM protocols. def __new__( cls, system=None, protocol=None, reference_system=None, exe=None, name="openmm", platform="CPU", work_dir=None, seed=None, property_map={}, **kwargs, ): from ._atm import OpenMMATM from ..Protocol._atm import _ATM # would like to use issubclass but _Protocol._ATM is not exposed if isinstance(protocol, _ATM): return super().__new__(OpenMMATM) else: return super().__new__(cls)
[docs] def __init__( self, system, protocol, reference_system=None, exe=None, name="openmm", platform="CPU", work_dir=None, seed=None, property_map={}, **kwargs, ): """ Constructor. Parameters ---------- system : :class:`System <BioSimSpace._SireWrappers.System>` The molecular system. protocol : :class:`Protocol <BioSimSpace.Protocol>` The protocol for the OpenMM process. reference_system : :class:`System <BioSimSpace._SireWrappers.System>` or None An optional system to use as a source of reference coordinates for position restraints. It is assumed that this system has the same topology as "system". If this is None, then "system" is used as a reference. exe : str The full path to the Python interpreter used to run OpenMM. name : str The name of the process. platform : str The platform for the simulation: "CPU", "CUDA", or "OPENCL". For CUDA use the CUDA_VISIBLE_DEVICES environment variable to set the GPUs on which to run, e.g. to run on two GPUs indexed 0 and 1 use: CUDA_VISIBLE_DEVICES=0,1. For OPENCL, instead use OPENCL_VISIBLE_DEVICES. work_dir : The working directory for the process. seed : int A random number seed. property_map : dict A dictionary that maps system "properties" to their user defined values. This allows the user to refer to properties with their own naming scheme, e.g. { "charge" : "my-charge" } kwargs : dict Additional keyword arguments. """ # Call the base class constructor. super().__init__( system, protocol, reference_system=reference_system, name=name, work_dir=work_dir, seed=seed, property_map=property_map, ) # Set the package name. self._package_name = "OPENMM" # This process can generate trajectory data. self._has_trajectory = True # If the path a Python interpreter wasn't specified, then use the bundled sire_python. if exe is None: bin_dir = _SireBase.getBinDir() # Generate the name of the sire_python interpreter. if _sys.platform == "win32": self._exe = _os.path.join(_os.path.normpath(bin_dir), "sire_python.exe") else: self._exe = _os.path.join(_os.path.normpath(bin_dir), "sire_python") else: # Make sure executable exists. if _os.path.isfile(exe): self._exe = exe else: raise IOError("OpenMM Python interpreter doesn't exist: '%s'" % exe) if not isinstance(platform, str): raise TypeError("'platform' must be of type 'str'.") else: # Strip all whitespace and convert to upper case. platform = platform.replace(" ", "").upper() # Check for platform support. if platform not in self._platforms: raise ValueError("Supported platforms are: %s" % self._platforms.keys()) else: self._platform = self._platforms[platform] # Initialise the stdout dictionary and title header. self._stdout_dict = _process._MultiDict() # Store the name of the OpenMM log file. self._log_file = _os.path.join(str(self._work_dir), f"{name}.log") # Initialise the log file separator. self._record_separator = None # Initialise a dictionary to map log file records to their column order. self._record_mapping = {} # The names of the input files. We choose to use AMBER files since they # are self-contained, but could equally work with GROMACS files. self._rst_file = _os.path.join(str(self._work_dir), f"{name}.rst7") self._top_file = _os.path.join(str(self._work_dir), f"{name}.prm7") self._ref_file = _os.path.join(str(self._work_dir), f"{name}_ref.rst7") # The name of the trajectory file. self._traj_file = _os.path.join(str(self._work_dir), f"{name}.dcd") # Set the path for the OpenMM Python script. (We use the concept of a # config file for consistency with other Process classes.) self._config_file = _os.path.join(str(self._work_dir), f"{name}_script.py") # Create the list of input files. self._input_files = [self._config_file, self._rst_file, self._top_file] # Add the reference file if there are position restraints. if self._protocol.getRestraint() is not None: self._input_files.append(self._ref_file) # Initialise the log file header. self._header = None # Now set up the working directory for the process. self._setup()
def __str__(self): """Return a human readable string representation of the object.""" return ( "<BioSimSpace.Process.%s: system=%s, protocol=%s, exe='%s', name='%s', platform='%s', work_dir='%s' seed=%s>" % ( self.__class__.__name__, str(self._system), self._protocol.__repr__(), self._exe + ("%s " % self._script if self._script else ""), self._name, self._platform, self._work_dir, self._seed, ) ) def __repr__(self): """Return a string showing how to instantiate the object.""" return ( "BioSimSpace.Process.%s(%s, %s, exe='%s', name='%s', platform='%s', work_dir='%s', seed=%s)" % ( self.__class__.__name__, str(self._system), self._protocol.__repr__(), self._exe + ("%s " % self._script if self._script else ""), self._name, self._platform, self._work_dir, self._seed, ) ) def _setup(self): """Setup the input files and working directory ready for simulation.""" # Create a copy of the system. system = self._system.copy() # Convert the water model topology so that it matches the AMBER naming convention. system._set_water_topology("AMBER", property_map=self._property_map) self._reference_system._set_water_topology( "AMBER", property_map=self._property_map ) # Check for perturbable molecules and convert to the chosen end state. system = self._checkPerturbable(system) reference_system = self._checkPerturbable(self._reference_system) # Create the input files... # RST file (coordinates). try: file = _os.path.splitext(self._rst_file)[0] _IO.saveMolecules(file, system, "rst7", property_map=self._property_map) except Exception as e: msg = "Failed to write system to 'RST7' format." if _isVerbose(): raise IOError(msg) from e else: raise IOError(msg) from None # Reference coordinate file for position restraints. if self._protocol.getRestraint() is not None: try: file = _os.path.splitext(self._ref_file)[0] _IO.saveMolecules( file, reference_system, "rst7", property_map=self._property_map, ) except Exception as e: msg = "Failed to write reference system to 'RST7' format." if _isVerbose(): raise IOError(msg) from e else: raise IOError(msg) from None # PRM file (topology). try: file = _os.path.splitext(self._top_file)[0] _IO.saveMolecules( file, system, "prm7", match_water=False, property_map=self._property_map ) except Exception as e: msg = "Failed to write system to 'PRM7' format." if _isVerbose(): raise IOError(msg) from e else: raise IOError(msg) from None # Skip if the user has passed a custom config. if isinstance(self._protocol, _Protocol.Custom): self.setConfig(self._protocol.getConfig()) else: self._generate_config() self.writeConfig(self._config_file) # Generate the dictionary of command-line arguments. self._generate_args() # Return the list of input files. return self._input_files def _generate_config(self): """Generate OpenMM Python script file strings.""" # Clear the existing configuration list. self._config = [] # Get the "space" property from the user mapping. prop = self._property_map.get("space", "space") # Check whether the system contains periodic box information. if prop in self._system._sire_object.propertyKeys(): try: # Make sure that we have a periodic box. The system will now have # a default cartesian space. box = self._system._sire_object.property(prop) has_box = box.isPeriodic() except: has_box = False else: _warnings.warn("No simulation box found. Assuming gas phase simulation.") has_box = False if isinstance(self._protocol, _Protocol.Minimisation): # Write the OpenMM import statements. self._add_config_imports() self._add_config_monkey_patches() # Load the input files. self.addToConfig("\n# Load the topology and coordinate files.") self.addToConfig( "\n# We use ParmEd due to issues with the built in AmberPrmtopFile for certain triclinic spaces." ) self.addToConfig( f"prm = parmed.load_file('{self._name}.prm7', '{self._name}.rst7')" ) # Don't use a cut-off if this is a vacuum simulation or if box information # is missing. self.addToConfig("\n# Initialise the molecular system.") if not has_box or not self._has_water: self.addToConfig("system = prm.createSystem(nonbondedMethod=NoCutoff,") else: self.addToConfig("system = prm.createSystem(nonbondedMethod=PME,") self.addToConfig(" nonbondedCutoff=1*nanometer,") self.addToConfig(" constraints=HBonds)") # Set the integrator. (Use zero-temperature as this is just a dummy step.) self.addToConfig("\n# Define the integrator.") self.addToConfig("integrator = LangevinMiddleIntegrator(0*kelvin,") self.addToConfig(" 1/picosecond,") self.addToConfig(" 0.002*picoseconds)") # Add the platform information. self._add_config_platform() # Add any position restraints. if self._protocol.getRestraint() is not None: self._add_config_restraints() # Set up the simulation object. self.addToConfig("\n# Initialise and configure the simulation object.") self.addToConfig("simulation = Simulation(prm.topology,") self.addToConfig(" system,") self.addToConfig(" integrator,") self.addToConfig(" platform,") self.addToConfig(" properties)") if self._protocol.getRestraint() is not None: self.addToConfig("simulation.context.setPositions(positions)") else: self.addToConfig("simulation.context.setPositions(prm.positions)") self.addToConfig("if prm.box_vectors is not None:") self.addToConfig( " box_vectors = reducePeriodicBoxVectors(prm.box_vectors)" ) self.addToConfig( " simulation.context.setPeriodicBoxVectors(*box_vectors)" ) self.addToConfig( f"simulation.minimizeEnergy(maxIterations={self._protocol.getSteps()})" ) # Add the reporters. self.addToConfig("\n# Add reporters.") self._add_config_reporters(state_interval=1, traj_interval=1) # Now run the simulation. self.addToConfig( "\n# Run a single simulation step to allow us to get the system and energy." ) self.addToConfig(f"simulation.step(1)") elif isinstance(self._protocol, _Protocol.Equilibration): # Write the OpenMM import statements and monkey-patches. self._add_config_imports() self._add_config_monkey_patches() # Load the input files. self.addToConfig("\n# Load the topology and coordinate files.") self.addToConfig( "\n# We use ParmEd due to issues with the built in AmberPrmtopFile for certain triclinic spaces." ) self.addToConfig( f"prm = parmed.load_file('{self._name}.prm7', '{self._name}.rst7')" ) # Don't use a cut-off if this is a vacuum simulation or if box information # is missing. is_periodic = True self.addToConfig("\n# Initialise the molecular system.") if not has_box or not self._has_water: is_periodic = False self.addToConfig("system = prm.createSystem(nonbondedMethod=NoCutoff,") else: self.addToConfig("system = prm.createSystem(nonbondedMethod=PME,") self.addToConfig(" nonbondedCutoff=1*nanometer,") self.addToConfig(" constraints=HBonds)") # Get the starting temperature and system pressure. temperature = self._protocol.getStartTemperature().kelvin().value() pressure = self._protocol.getPressure() # Add a Monte Carlo barostat if the simulation is at constant pressure. is_const_pressure = False if pressure is not None: # Cannot use a barostat with a non-periodic system. if not is_periodic: _warnings.warn( "Cannot use a barostat for a vacuum or non-periodic simulation" ) else: is_const_pressure = True # Convert to bar and get the value. pressure = pressure.bar().value() # Create the barostat and add its force to the system. self.addToConfig("\n# Add a barostat to run at constant pressure.") self.addToConfig( f"barostat = MonteCarloBarostat({pressure}*bar, {temperature}*kelvin)" ) if self._is_seeded: self.addToConfig(f"barostat.setRandomNumberSeed({self._seed})") self.addToConfig("system.addForce(barostat)") # Add any position restraints. if self._protocol.getRestraint() is not None: self._add_config_restraints() # Get the integration time step from the protocol. timestep = self._protocol.getTimeStep().picoseconds().value() # Set the integrator. self.addToConfig("\n# Define the integrator.") self.addToConfig( f"integrator = LangevinMiddleIntegrator({temperature}*kelvin," ) friction = ( 1 / self._protocol.getThermostatTimeConstant().picoseconds().value() ) self.addToConfig( f" {friction:.5f}/picosecond," ) self.addToConfig(f" {timestep}*picoseconds)") if self._is_seeded: self.addToConfig(f"integrator.setRandomNumberSeed({self._seed})") # Add the platform information. self._add_config_platform() # Set up the simulation object. self.addToConfig("\n# Initialise and configure the simulation object.") self.addToConfig("simulation = Simulation(prm.topology,") self.addToConfig(" system,") self.addToConfig(" integrator,") self.addToConfig(" platform,") self.addToConfig(" properties)") if self._protocol.getRestraint() is not None: self.addToConfig("simulation.context.setPositions(positions)") else: self.addToConfig("simulation.context.setPositions(prm.positions)") self.addToConfig("if prm.box_vectors is not None:") self.addToConfig( " box_vectors = reducePeriodicBoxVectors(prm.box_vectors)" ) self.addToConfig( " simulation.context.setPeriodicBoxVectors(*box_vectors)" ) # Set initial velocities from temperature distribution. self.addToConfig("\n# Setting initial system velocities.") self.addToConfig( f"simulation.context.setVelocitiesToTemperature({temperature})" ) # Work out the number of integration steps. steps = _math.ceil( self._protocol.getRunTime() / self._protocol.getTimeStep() ) # Get the report and restart intervals. report_interval = self._protocol.getReportInterval() restart_interval = self._protocol.getRestartInterval() # Cap the intervals at the total number of steps. if report_interval > steps: report_interval = steps if restart_interval > steps: restart_interval = steps # Add the reporters. self.addToConfig("\n# Add reporters.") self._add_config_reporters( state_interval=report_interval, traj_interval=restart_interval, is_restart=False, ) # Now run the simulation. self.addToConfig("\n# Run the simulation.") # Constant temperature equilibration. if self._protocol.isConstantTemp(): self.addToConfig(f"simulation.step({steps})") # Heating / cooling cycle. else: # Adjust temperature every 100 cycles, assuming that there at # least that many cycles. if steps > 100: # Work out the number of temperature cycles. temp_cycles = _math.ceil(steps / 100) # Work out the temperature change per cycle. delta_temp = ( self._protocol.getEndTemperature().kelvin().value() - self._protocol.getStartTemperature().kelvin().value() ) / temp_cycles self.addToConfig(f"start_temperature = {temperature}") self.addToConfig(f"for x in range(0, {temp_cycles}):") self.addToConfig( f" temperature = {temperature} + x*{delta_temp}" ) self.addToConfig( f" integrator.setTemperature(temperature*kelvin)" ) if is_const_pressure: self.addToConfig( f" barostat.setDefaultTemperature(temperature*kelvin)" ) self.addToConfig(" simulation.step(100)") else: # Work out the temperature change per step. delta_temp = ( self._protocol.getEndTemperature().kelvin().value() - self._protocol.getStartTemperature().kelvin().value() ) / steps self.addToConfig(f"start_temperature = {temperature}") self.addToConfig(f"for x in range(0, {steps}):") self.addToConfig( f" temperature = {temperature} + x*{delta_temp}" ) self.addToConfig( f" integrator.setTemperature(temperature*kelvin)" ) if is_const_pressure: self.addToConfig( f" barostat.setDefaultTemperature(temperature*kelvin)" ) self.addToConfig(" simulation.step(1)") elif isinstance(self._protocol, _Protocol.Production): # Write the OpenMM import statements. self._add_config_imports() self._add_config_monkey_patches() # Production specific import. self.addToConfig("import os") # Load the input files. self.addToConfig("\n# Load the topology and coordinate files.") self.addToConfig( "\n# We use ParmEd due to issues with the built in AmberPrmtopFile for certain triclinic spaces." ) self.addToConfig( f"prm = parmed.load_file('{self._name}.prm7', '{self._name}.rst7')" ) # Don't use a cut-off if this is a vacuum simulation or if box information # is missing. is_periodic = True self.addToConfig("\n# Initialise the molecular system.") if not has_box or not self._has_water: is_periodic = False self.addToConfig("system = prm.createSystem(nonbondedMethod=NoCutoff,") else: self.addToConfig("system = prm.createSystem(nonbondedMethod=PME,") self.addToConfig(" nonbondedCutoff=1*nanometer,") self.addToConfig(" constraints=HBonds)") # Get the starting temperature and system pressure. temperature = self._protocol.getTemperature().kelvin().value() pressure = self._protocol.getPressure() # Add a Monte Carlo barostat if the simulation is at constant pressure. is_const_pressure = False if pressure is not None: # Cannot use a barostat with a non-periodic system. if not is_periodic: _warnings.warn( "Cannot use a barostat for a vacuum or non-periodic simulation" ) else: is_const_pressure = True # Convert to bar and get the value. pressure = pressure.bar().value() # Create the barostat and add its force to the system. self.addToConfig("\n# Add a barostat to run at constant pressure.") self.addToConfig( f"barostat = MonteCarloBarostat({pressure}*bar, {temperature}*kelvin)" ) if self._is_seeded: self.addToConfig(f"barostat.setRandomNumberSeed({self._seed})") self.addToConfig("system.addForce(barostat)") # Add any position restraints. if self._protocol.getRestraint() is not None: self._add_config_restraints() # Get the integration time step from the protocol. timestep = self._protocol.getTimeStep().picoseconds().value() # Set the integrator. self.addToConfig("\n# Define the integrator.") self.addToConfig( f"integrator = LangevinMiddleIntegrator({temperature}*kelvin," ) friction = ( 1 / self._protocol.getThermostatTimeConstant().picoseconds().value() ) self.addToConfig( f" {friction:.5f}/picosecond," ) self.addToConfig(f" {timestep}*picoseconds)") if self._is_seeded: self.addToConfig(f"integrator.setRandomNumberSeed({self._seed})") # Add the platform information. self._add_config_platform() # Set up the simulation object. self.addToConfig("\n# Initialise and configure the simulation object.") self.addToConfig("simulation = Simulation(prm.topology,") self.addToConfig(" system,") self.addToConfig(" integrator,") self.addToConfig(" platform,") self.addToConfig(" properties)") if self._protocol.getRestraint() is not None: self.addToConfig("simulation.context.setPositions(positions)") else: self.addToConfig("simulation.context.setPositions(prm.positions)") self.addToConfig("if prm.box_vectors is not None:") self.addToConfig( " box_vectors = reducePeriodicBoxVectors(prm.box_vectors)" ) self.addToConfig( " simulation.context.setPeriodicBoxVectors(*box_vectors)" ) # Set initial velocities from temperature distribution. self.addToConfig("\n# Setting initial system velocities.") self.addToConfig( f"simulation.context.setVelocitiesToTemperature({temperature})" ) # Check for a restart file and load the simulation state. is_restart, step = self._add_config_restart() # Work out the number of integration steps. total_steps = _math.ceil( self._protocol.getRunTime() / self._protocol.getTimeStep() ) # Subtract the current number of steps. steps = total_steps - step # Exit if the simulation has already finished. if steps <= 0: print("The simulation has already finished!") return # Inform user that a restart was loaded. self.addToConfig("\n# Print restart information.") self.addToConfig("if is_restart:") self.addToConfig(f" steps = {total_steps}") self.addToConfig(" percent_complete = 100 * (step / steps)") self.addToConfig(" print('Loaded state from an existing simulation.')") self.addToConfig( " print(f'Simulation is {percent_complete}% complete.')" ) # Get the report and restart intervals. report_interval = self._protocol.getReportInterval() restart_interval = self._protocol.getRestartInterval() # Cap the intervals at the total number of steps. if report_interval > steps: report_interval = steps if restart_interval > steps: restart_interval = steps # Add the reporters. self.addToConfig("\n# Add reporters.") self._add_config_reporters( state_interval=report_interval, traj_interval=restart_interval, is_restart=is_restart, ) # Work out the total simulation time in picoseconds. run_time = steps * timestep # Work out the number of cycles in 100 picosecond intervals. cycles = _math.ceil(run_time / 100) # Work out the number of steps per cycle. steps_per_cycle = int(steps / cycles) # Now run the simulation. self.addToConfig("\n# Run the simulation in 100 picosecond cycles.") self.addToConfig(f"for x in range(0, {cycles}):") self.addToConfig(f" simulation.step({steps_per_cycle})") self.addToConfig(f" simulation.saveState('{self._name}.xml')") elif isinstance(self._protocol, _Protocol.Metadynamics): colvar = self._protocol.getCollectiveVariable() if len(colvar) != 1 or not isinstance( colvar[0], _CollectiveVariable.Funnel ): raise _IncompatibleError( "We currently only support '%s' collective variables for '%s' protocols" % ( _CollectiveVariable.Funnel.__name__, self._protocol.__class__.__name__, ) ) # Create the path to the patched OpenMM metadynamics module. aux_file = "metadynamics.py" path = ( _os.path.dirname(_CollectiveVariable.__file__).replace( "CollectiveVariable", "_aux" ) + "/" + aux_file ) # Copy the file into the working directory. _shutil.copyfile(path, _os.path.join(str(self._work_dir), aux_file)) # The following OpenMM native implementation of the funnel metadynamics protocol # is adapted from funnel_maker.py by Dominykas Lukauskis. # Extract the only collective variable. colvar = colvar[0] # Write the OpenMM import statements. self._add_config_imports() self.addToConfig( "from metadynamics import *" ) # Use local patched metadynamics module. self.addToConfig("from glob import glob") self.addToConfig("import math") self.addToConfig("import os") self.addToConfig("import shutil") # Load the input files. self.addToConfig("\n# Load the topology and coordinate files.") self.addToConfig( "\n# We use ParmEd due to issues with the built in AmberPrmtopFile for certain triclinic spaces." ) self.addToConfig( f"prm = parmed.load_file('{self._name}.prm7', '{self._name}.rst7')" ) # Don't use a cut-off if this is a vacuum simulation or if box information # is missing. is_periodic = True self.addToConfig("\n# Initialise the molecular system.") if not has_box or not self._has_water: is_periodic = False self.addToConfig("system = prm.createSystem(nonbondedMethod=NoCutoff,") else: self.addToConfig("system = prm.createSystem(nonbondedMethod=PME,") self.addToConfig(" nonbondedCutoff=1*nanometer,") self.addToConfig(" constraints=HBonds)") # Get the starting temperature and system pressure. temperature = self._protocol.getTemperature().kelvin().value() pressure = self._protocol.getPressure() # Add a Monte Carlo barostat if the simulation is at constant pressure. is_const_pressure = False if pressure is not None: # Cannot use a barostat with a non-periodic system. if not is_periodic: _warnings.warn( "Cannot use a barostat for a vacuum or non-periodic simulation" ) else: is_const_pressure = True # Convert to bar and get the value. pressure = pressure.bar().value() # Create the barostat and add its force to the system. self.addToConfig("\n# Add a barostat to run at constant pressure.") self.addToConfig( f"barostat = MonteCarloBarostat({pressure}*bar, {temperature}*kelvin)" ) if self._is_seeded: self.addToConfig(f"barostat.setRandomNumberSeed({self._seed})") self.addToConfig("system.addForce(barostat)") # Store the number of atoms in each molecule. atom_nums = [] for mol in self._system: atom_nums.append(mol.nAtoms()) # Sort the molecule indices by the number of atoms they contain. sorted_nums = sorted((value, idx) for idx, value in enumerate(atom_nums)) # Set the ligand to the index of the second largest molecule. ligand = sorted_nums[-2][1] # Work out the start/end indices of the ligand within the system. idx_start = self._system.getIndex(self._system[ligand].getAtoms()[0]) idx_end = self._system.getIndex(self._system[ligand + 1].getAtoms()[0]) # Add the funnel variables. p1_string = ",".join([str(x) for x in colvar.getAtoms0()]) p2_string = ",".join([str(x) for x in colvar.getAtoms1()]) self.addToConfig("\n# Set funnel variables.") self.addToConfig(f"p1 = [{p1_string}]") self.addToConfig(f"p2 = [{p2_string}]") self.addToConfig(f"lig = [x for x in range({idx_start}, {idx_end})]") sigma_proj = colvar.getHillWidth()[0].nanometers().value() self.addToConfig("\n# Create the bias variable for the funnel projection.") self.addToConfig( "projection = CustomCentroidBondForce(3, 'distance(g1,g2)*cos(angle(g1,g2,g3))')" ) self.addToConfig("projection.addGroup(lig)") self.addToConfig("projection.addGroup(p1)") self.addToConfig("projection.addGroup(p2)") self.addToConfig("projection.addBond([0,1,2])") self.addToConfig("projection.setUsesPeriodicBoundaryConditions(True)") self.addToConfig(f"sigma_proj = {sigma_proj}") if colvar.getLowerBound() is None and colvar.getUpperBound() is None: # Sane defaults if no bounds are set. lower_wall = 0.5 upper_wall = 4.5 else: lower_wall = colvar.getLowerBound().getValue().nanometers().value() upper_wall = colvar.getUpperBound().getValue().nanometers().value() self.addToConfig( f"proj = BiasVariable(projection, {lower_wall-0.2}, {upper_wall+0.2}, {sigma_proj}, False, gridWidth=200)" ) sigma_ext = colvar.getHillWidth()[1].nanometers().value() self.addToConfig("\n# Create the bias variable for the funnel extent.") self.addToConfig( "extent = CustomCentroidBondForce(3, 'distance(g1,g2)*sin(angle(g1,g2,g3))')" ) self.addToConfig("extent.addGroup(lig)") self.addToConfig("extent.addGroup(p1)") self.addToConfig("extent.addGroup(p2)") self.addToConfig("extent.addBond([0,1,2])") self.addToConfig("extent.setUsesPeriodicBoundaryConditions(True)") extent_max = ( colvar.getWidth().nanometers().value() + colvar.getBuffer().nanometers().value() + 0.2 ) self.addToConfig(f"sigma_ext = {sigma_ext}") self.addToConfig( f"ext = BiasVariable(extent, 0.0, {extent_max}, {sigma_ext}, False, gridWidth=200)" ) # Add restraints. self.addToConfig("\n# Add restraints.") self.addToConfig("k1 = 10000*kilojoules_per_mole") self.addToConfig("k2 = 1000*kilojoules_per_mole") self.addToConfig(f"lower_wall = {lower_wall}*nanometer") self.addToConfig(f"upper_wall = {upper_wall}*nanometer") self.addToConfig("\n# Upper wall.") self.addToConfig( "upper_wall_rest = CustomCentroidBondForce(3, '(k1/2)*max(distance(g1,g2)*cos(angle(g1,g2,g3)) - upper_wall, 0)^2')" ) self.addToConfig("upper_wall_rest.addGroup(lig)") self.addToConfig("upper_wall_rest.addGroup(p1)") self.addToConfig("upper_wall_rest.addGroup(p2)") self.addToConfig("upper_wall_rest.addBond([0,1,2])") self.addToConfig("upper_wall_rest.addGlobalParameter('k1', k1)") self.addToConfig( "upper_wall_rest.addGlobalParameter('upper_wall', upper_wall)" ) self.addToConfig("upper_wall_rest.setUsesPeriodicBoundaryConditions(True)") self.addToConfig("system.addForce(upper_wall_rest)") self.addToConfig("\n# Sides of the funnel.") self.addToConfig( f"wall_width = {colvar.getWidth().nanometers().value()}*nanometer" ) self.addToConfig( f"wall_buffer = {colvar.getBuffer().nanometers().value()}*nanometer" ) self.addToConfig(f"beta_cent = {colvar.getSteepness()}") self.addToConfig( f"s_cent = {colvar.getInflection().nanometers().value()}*nanometer" ) self.addToConfig( "dist_restraint = CustomCentroidBondForce(3, '(k2/2)*max(distance(g1,g2)*sin(angle(g1,g2,g3)) - (a/(1+exp(b*(distance(g1,g2)*cos(angle(g1,g2,g3))-c)))+d), 0)^2')" ) self.addToConfig("dist_restraint.addGroup(lig)") self.addToConfig("dist_restraint.addGroup(p1)") self.addToConfig("dist_restraint.addGroup(p2)") self.addToConfig("dist_restraint.addBond([0,1,2])") self.addToConfig("dist_restraint.addGlobalParameter('k2', k2)") self.addToConfig("dist_restraint.addGlobalParameter('a', wall_width)") self.addToConfig("dist_restraint.addGlobalParameter('b', beta_cent)") self.addToConfig("dist_restraint.addGlobalParameter('c', s_cent)") self.addToConfig("dist_restraint.addGlobalParameter('d', wall_buffer)") self.addToConfig("dist_restraint.setUsesPeriodicBoundaryConditions(True)") self.addToConfig("system.addForce(dist_restraint)") self.addToConfig("\n# Lower wall.") self.addToConfig( "lower_wall_rest = CustomCentroidBondForce(3, '(k1/2)*min(distance(g1,g2)*cos(angle(g1,g2,g3)) - lower_wall, 0)^2')" ) self.addToConfig("lower_wall_rest.addGroup(lig)") self.addToConfig("lower_wall_rest.addGroup(p1)") self.addToConfig("lower_wall_rest.addGroup(p2)") self.addToConfig("lower_wall_rest.addBond([0,1,2])") self.addToConfig("lower_wall_rest.addGlobalParameter('k1', k1)") self.addToConfig( "lower_wall_rest.addGlobalParameter('lower_wall', lower_wall)" ) self.addToConfig("lower_wall_rest.setUsesPeriodicBoundaryConditions(True)") self.addToConfig("system.addForce(lower_wall_rest)") # Get the number of steps to date. step = 0 if _os.path.isfile(_os.path.join(str(self._work_dir), f"{self._name}.xml")): if _os.path.isfile( _os.path.join(str(self._work_dir), f"{self._name}.log") ): with open( _os.path.join(str(self._work_dir), f"{self._name}.log"), "r" ) as f: lines = f.readlines() last_line = lines[-1].split() try: step = int(last_line[0]) except: raise _IncompatibleError( "Failed to read current integration " f"step from '{self._name}.log'" ) else: raise IOError(f"Missing log file: '{self._name}.log'") # Get the report and restart intervals. report_interval = self._protocol.getReportInterval() restart_interval = self._protocol.getRestartInterval() # Get the hill deposition frequency. hill_freq = self._protocol.getHillFrequency() # Work out the number of integration steps. total_steps = _math.ceil( self._protocol.getRunTime() / self._protocol.getTimeStep() ) # Work out the number of cycles. # (Record COLVAR and HILLS at hill deposition frequency.) total_cycles = _math.ceil(total_steps / hill_freq) # Subtract the current number of steps. remaining_steps = total_steps - step # Exit if the simulation has already finished. if remaining_steps <= 0: print("The simulation has already finished!") return # Cap the intervals at the total number of steps. if report_interval > remaining_steps: report_interval = remaining_steps if restart_interval > remaining_steps: restart_interval = remaining_steps self.addToConfig("\n# Initialise the metadynamics object.") if self._protocol.getBiasFactor() is None: bias = 1.0000001 else: bias = self._protocol.getBiasFactor() self.addToConfig(f"bias = {bias}") height = self._protocol.getHillHeight().kj_per_mol().value() self.addToConfig( f"meta = Metadynamics(system, [proj, ext], {temperature}*kelvin, bias, {height}*kilojoules_per_mole, {hill_freq}, biasDir = '.', saveFrequency = {report_interval})" ) # Get the integration time step from the protocol. timestep = self._protocol.getTimeStep().picoseconds().value() # Set the integrator. self.addToConfig("\n# Define the integrator.") self.addToConfig( f"integrator = LangevinMiddleIntegrator({temperature}*kelvin," ) self.addToConfig(" 1/picosecond,") self.addToConfig(f" {timestep}*picoseconds)") if self._is_seeded: self.addToConfig(f"integrator.setRandomNumberSeed({self._seed})") # Add the platform information. self._add_config_platform() # Set up the simulation object. self.addToConfig("\n# Initialise and configure the simulation object.") self.addToConfig("simulation = Simulation(prm.topology,") self.addToConfig(" system,") self.addToConfig(" integrator,") self.addToConfig(" platform,") self.addToConfig(" properties)") self.addToConfig("simulation.context.setPositions(prm.positions)") self.addToConfig("if prm.box_vectors is not None:") self.addToConfig( " box_vectors = reducePeriodicBoxVectors(prm.box_vectors)" ) self.addToConfig( " simulation.context.setPeriodicBoxVectors(*box_vectors)" ) # Set initial velocities from temperature distribution. self.addToConfig("\n# Setting initial system velocities.") self.addToConfig( f"simulation.context.setVelocitiesToTemperature({temperature})" ) # Check for a restart file and load the simulation state. is_restart, step = self._add_config_restart() # Add the reporters. self.addToConfig("\n# Add reporters.") self._add_config_reporters( state_interval=report_interval, traj_interval=restart_interval, is_restart=is_restart, ) # Create the HILLS file. self.addToConfig("\n# Create PLUMED compatible HILLS file.") self.addToConfig("if is_restart:") self.addToConfig(" file = open('HILLS','a')") self.addToConfig("else:") self.addToConfig(" file = open('HILLS','w')") self.addToConfig( " file.write('#! FIELDS time pp.proj pp.ext sigma_pp.proj sigma_pp.ext height biasf\\n')" ) self.addToConfig(" file.write('#! SET multivariate false\\n')") self.addToConfig(" file.write('#! SET kerneltype gaussian\\n')") # Get the initial collective variables. self.addToConfig("\n# Initialise the collective variable array.") self.addToConfig( "current_cvs = np.array(list(meta.getCollectiveVariables(simulation)) + [meta.getHillHeight(simulation)])" ) # Write the initial record. self.addToConfig("\n# Write the initial collective variable record.") self.addToConfig("if is_restart:") self.addToConfig(" if os.path.isfile('COLVAR.npy'):") self.addToConfig(" colvar_array = np.load('COLVAR.npy')") self.addToConfig( " colvar_array = np.append(colvar_array, [current_cvs], axis=0)" ) self.addToConfig(" else:") self.addToConfig(" raise IOError('Missing COLVAR file: COLVAR.npy')") self.addToConfig("else:") self.addToConfig(" colvar_array = np.array([current_cvs])") self.addToConfig(" line = colvar_array[0]") self.addToConfig(" time = 0") self.addToConfig( " write_line = f'{time:15} {line[0]:20.16f} {line[1]:20.16f} {sigma_proj} {sigma_ext} {line[2]:20.16f} {bias}\\n'" ) self.addToConfig(" file.write(write_line)") # Run the metadynamics simulation. self.addToConfig("\n# Run the simulation.") self.addToConfig(f"total_steps = {total_steps}") self.addToConfig(f"total_cycles = {total_cycles}") self.addToConfig(f"remaining_steps = {remaining_steps}") self.addToConfig("steps_per_cycle = math.ceil(total_steps / total_cycles)") self.addToConfig( "remaining_cycles = math.ceil(remaining_steps / steps_per_cycle)" ) self.addToConfig(f"start_cycles = total_cycles - remaining_cycles") self.addToConfig("checkpoint = 100") self.addToConfig("if is_restart:") self.addToConfig(" fraction_complete = step / total_steps") self.addToConfig(" percent_complete = 100 * fraction_complete") self.addToConfig(" print('Loaded state from an existing simulation.')") self.addToConfig( " print(f'Simulation is {percent_complete}% complete.')" ) self.addToConfig(" if fraction_complete > 1.0:") self.addToConfig( " start_cycles = math.ceil(fraction_complete * total_cycles)" ) self.addToConfig(" total_cycles = start_cycles + total_cycles") self.addToConfig("for x in range(start_cycles, total_cycles):") self.addToConfig(" meta.step(simulation, steps_per_cycle)") self.addToConfig( " current_cvs = np.array(list(meta.getCollectiveVariables(simulation)) + [meta.getHillHeight(simulation)])" ) self.addToConfig( " colvar_array = np.append(colvar_array, [current_cvs], axis=0)" ) self.addToConfig(" np.save('COLVAR.npy', colvar_array)") self.addToConfig(" line = colvar_array[x+1]") self.addToConfig(f" time = int((x+1) * {timestep}*steps_per_cycle)") self.addToConfig( " write_line = f'{time:15} {line[0]:20.16f} {line[1]:20.16f} {sigma_proj} {sigma_ext} {line[2]:20.16f} {bias}\\n'" ) self.addToConfig(" file.write(write_line)") self.addToConfig(" # Save the simulation state every 100 picoseconds.") self.addToConfig(" if time >= checkpoint:") self.addToConfig(f" simulation.saveState('{self._name}.xml')") self.addToConfig(" while time >= checkpoint:") self.addToConfig(" checkpoint += 100") # Create a dummy PLUMED input file so that we can bind PLUMED # analysis functions to this process. self._plumed = _Plumed(str(self._work_dir)) plumed_config, auxillary_files = self._plumed.createConfig( self._system, self._protocol, self._property_map ) # Expose the PLUMED specific member functions. setattr(self, "getFreeEnergy", self._getFreeEnergy) setattr(self, "getCollectiveVariable", self._getCollectiveVariable) setattr(self, "sampleConfigurations", self._sampleConfigurations) setattr(self, "getTime", self._getTime) else: raise _IncompatibleError( "Unsupported protocol: '%s'" % self._protocol.__class__.__name__ ) # Flag that this isn't a custom protocol. self._protocol._setCustomised(False)
[docs] def getConfig(self): """ Get the list of strings defining the OpenMM Python script. Returns ------- config : [str] The list of configuration strings. """ return super().getConfig()
def setConfig(self, config): """ Set the list of strings defining the OpenMM Python script. Parameters ---------- config : str, [str] The list of configuration strings, or a path to a configuration file. """ return super().setConfig()
[docs] def addToConfig(self, config): """ Add a string to the OpenMM Python script configuration. Parameters ---------- config : str, [str] A configuration string, a list of configuration strings, or a path to a configuration file. """ # Call the base class method. super().addToConfig(config)
[docs] def resetConfig(self): """Reset the OpenMM Python script configuration.""" self._generate_config()
[docs] def setConfig(self, config): """ Set the list of configuration file strings. Parameters ---------- config : str, [str] The list of configuration strings, or a path to a configuration file. """ # Call the base class method. super().setConfig(config)
[docs] def start(self): """ Start the OpenMM process. Returns ------- process : :class:`Process.OpenMM <BioSimSpace.Process.OpenMM>` A handle to the OpenMM process. """ # The process is currently queued. if self.isQueued(): return # Process is already running. if self._process is not None: if self._process.isRunning(): return # Clear any existing output. self._clear_output() # Run the process in the working directory. with _Utils.cd(self._work_dir): # Create the arguments string list. # The name of the Python script (config file) is the first argument. args = ["%s" % self._config_file] args.extend(self.getArgStringList()) # Write the command-line process to a README.txt file. with open("README.txt", "w") as f: # Set the command-line string. self._command = ( "%s %s " % (self._exe, self._config_file) + self.getArgString() ) # Write the command to file. f.write("# OpenMM was run with the following command:\n") f.write("%s\n" % self._command) # Start the timer. self._timer = _timeit.default_timer() # Start the simulation. self._process = _SireBase.Process.run( self._exe, args, "%s.out" % self._name, "%s.err" % self._name ) return self
[docs] def getSystem(self, block="AUTO"): """ Get the latest molecular system. Parameters ---------- block : bool Whether to block until the process has finished running. Returns ------- system : :class:`System <BioSimSpace._SireWrappers.System>` The latest molecular system. """ # Wait for the process to finish. if block is True: self.wait() elif block == "AUTO" and self._is_blocked: self.wait() # Try to get the most recent trajectory frame. try: # Handle minimisation protocols separately. if isinstance( self._protocol, (_Protocol.Minimisation, _Protocol.ATMMinimisation) ): # Do we need to get coordinates for the lambda=1 state. if "is_lambda1" in self._property_map: is_lambda1 = True else: is_lambda1 = False traj = self.getTrajectory() # If there is no trajectory, simply return None. if traj is None: return None # Get the last frame. new_system = traj.getFrames(-1)[0] # Make a copy of the existing molecular system. old_system = self._system.copy() # Update the coordinates and velocities and return a mapping between # the molecule indices in the two systems. sire_system, mapping = _SireIO.updateCoordinatesAndVelocities( old_system._sire_object, new_system._sire_object, self._mapping, is_lambda1, self._property_map, self._property_map, ) # Update the underlying Sire object. old_system._sire_object = sire_system # Store the mapping between the MolIdx in both systems so we don't # need to recompute it next time. self._mapping = mapping # Update the box information in the original system. if "space" in new_system._sire_object.propertyKeys(): box = new_system._sire_object.property("space") if box.isPeriodic(): old_system._sire_object.setProperty( self._property_map.get("space", "space"), box ) return old_system else: # Work out the total number of trajectory frames. num_frames = int( (self._protocol.getRunTime() / self._protocol.getTimeStep()) / self._protocol.getRestartInterval() ) # Work out the fraction of the simulation that has been completed. frac_complete = self._protocol.getRunTime() / self.getTime() # Make sure the fraction doesn't exceed one. OpenMM can report # time values that are larger than the number of integration steps # multiplied by the time step. if frac_complete > 1: frac_complete = 1 # Work out the trajectory frame index, rounding down. # Remember that frames in MDTraj are zero indexed, like Python. index = int(frac_complete * num_frames) if index > 0: index -= 1 # Return the most recent frame. return self.getFrame(index) except: return None
[docs] def getCurrentSystem(self): """ Get the latest molecular system. Returns ------- system : :class:`System <BioSimSpace._SireWrappers.System>` The latest molecular system. """ return self.getSystem(block=False)
[docs] def getTrajectory(self, backend="AUTO", block="AUTO"): """ Return a trajectory object. Parameters ---------- backend : str The backend to use for trajectory parsing. To see supported backends, run BioSimSpace.Trajectory.backends(). Using "AUTO" will try each in sequence. block : bool Whether to block until the process has finished running. Returns ------- trajectory : :class:`System <BioSimSpace.Trajectory.Trajectory>` The latest trajectory object. """ if not isinstance(backend, str): raise TypeError("'backend' must be of type 'str'") if not isinstance(block, (bool, str)): raise TypeError("'block' must be of type 'bool' or 'str'") # Wait for the process to finish. if block is True: self.wait() elif block == "AUTO" and self._is_blocked: self.wait() # Warn the user if the process has exited with an error. if self.isError(): _warnings.warn("The process exited with an error!") if not _os.path.isfile(self._traj_file): return None else: return _Trajectory.Trajectory(process=self, backend=backend)
[docs] def getFrame(self, index): """ Return a specific trajectory frame. Parameters ---------- index : int The index of the frame. Returns ------- frame : :class:`System <BioSimSpace._SireWrappers.System>` The System object of the corresponding frame. """ 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 ) if index < 0 or index > max_index: raise ValueError(f"'index' must be in range [0, {max_index}].") # Warn the user if the process has exited with an error. if self.isError(): _warnings.warn("The process exited with an error!") try: # Do we need to get coordinates for the lambda=1 state. if "is_lambda1" in self._property_map: is_lambda1 = True else: is_lambda1 = False # Create a copy of the existing system object. old_system = self._system.copy() # Get the latest trajectory frame. new_system = _Trajectory.getFrame(self._traj_file, self._top_file, index) # Update the coordinates and velocities and return a mapping between # the molecule indices in the two systems. sire_system, mapping = _SireIO.updateCoordinatesAndVelocities( old_system._sire_object, new_system._sire_object, self._mapping, is_lambda1, self._property_map, self._property_map, ) # Update the underlying Sire object. old_system._sire_object = sire_system # Store the mapping between the MolIdx in both systems so we don't # need to recompute it next time. self._mapping = mapping # Update the box information in the original system. if "space" in new_system._sire_object.propertyKeys(): box = new_system._sire_object.property("space") old_system._sire_object.setProperty( self._property_map.get("space", "space"), box ) return old_system except: return None
[docs] def getRecord(self, record, time_series=False, unit=None, block="AUTO"): """ Get a record from the stdout dictionary. Parameters ---------- record : str The record key. time_series : bool Whether to return a list of time series records. unit : :class:`Unit <BioSimSpace.Units>` The unit to convert the record to. block : bool Whether to block until the process has finished running. Returns ------- record : :class:`Type <BioSimSpace.Types>` The matching record. """ # Wait for the process to finish. if block is True: self.wait() elif block == "AUTO" and self._is_blocked: self.wait() # Warn the user if the process has exited with an error. if self.isError(): _warnings.warn("The process exited with an error!") self._update_stdout_dict() return self._get_stdout_record(record, time_series, unit)
[docs] def getCurrentRecord(self, record, time_series=False, unit=None): """ Get a current record from the stdout dictionary. Parameters ---------- record : str The record key. time_series : bool Whether to return a list of time series records. unit : :class:`Unit <BioSimSpace.Units>` The unit to convert the record to. Returns ------- record : :class:`Type <BioSimSpace.Types>` The matching record. """ # Warn the user if the process has exited with an error. if self.isError(): _warnings.warn("The process exited with an error!") self._update_stdout_dict() return self._get_stdout_record(record, time_series, unit)
[docs] def getRecords(self, block="AUTO"): """ Return the dictionary of stdout time-series records. Parameters ---------- block : bool Whether to block until the process has finished running. Returns ------- records : :class:`MultiDict <BioSimSpace.Process._process._MultiDict>` The dictionary of time-series records. """ # Wait for the process to finish. if block is True: self.wait() elif block == "AUTO" and self._is_blocked: self.wait() # Warn the user if the process has exited with an error. if self.isError(): _warnings.warn("The process exited with an error!") self._update_stdout_dict() return self._stdout_dict.copy()
[docs] def getCurrentRecords(self): """ Return the current dictionary of stdout time-series records. Parameters ---------- block : bool Whether to block until the process has finished running. Returns ------- records : :class:`MultiDict <BioSimSpace.Process._process._MultiDict>` The dictionary of time-series records. """ return self.getRecords(block=False)
[docs] def getTime(self, time_series=False, block="AUTO"): """ Get the simulation time. Parameters ---------- time_series : bool Whether to return a list of time series records. block : bool Whether to block until the process has finished running. Returns ------- time : :class:`Time <BioSimSpace.Types.Time>` The current simulation time in nanoseconds. """ return self.getRecord("TIME(PS)", time_series, _Units.Time.picosecond, block)
[docs] def getCurrentTime(self, time_series=False): """ Get the current simulation time. Parameters ---------- time_series : bool Whether to return a list of time series records. block : bool Whether to block until the process has finished running. Returns ------- time : :class:`Time <BioSimSpace.Types.Time>` The current simulation time in nanoseconds. """ return self.getTime(time_series, block=False)
[docs] def getStep(self, time_series=False, block="AUTO"): """ Get the number of integration steps. Parameters ---------- time_series : bool Whether to return a list of time series records. block : bool Whether to block until the process has finished running. Returns ------- step : int The current number of integration steps. """ return self.getRecord("STEP", time_series, None, block)
[docs] def getCurrentStep(self, time_series=False): """ Get the current number of integration steps. Parameters ---------- time_series : bool Whether to return a list of time series records. Returns ------- step : int The current number of integration steps. """ return self.getStep(time_series, block=False)
[docs] def getPotentialEnergy(self, time_series=False, block="AUTO"): """ Get the potential energy. Parameters ---------- time_series : bool Whether to return a list of time series records. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The potential energy. """ return self.getRecord( "POTENTIALENERGY(KJ/MOLE)", time_series, _Units.Energy.kj_per_mol, block )
[docs] def getCurrentPotentialEnergy(self, time_series=False): """ Get the current potential energy. Parameters ---------- time_series : bool Whether to return a list of time series records. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The potential energy. """ return self.getPotentialEnergy(time_series, block=False)
[docs] def getKineticEnergy(self, time_series=False, block="AUTO"): """ Get the kinetic energy. Parameters ---------- time_series : bool Whether to return a list of time series records. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The kinetic energy. """ return self.getRecord( "KINETICENERGY(KJ/MOLE)", time_series, _Units.Energy.kj_per_mol, block )
[docs] def getCurrentKineticEnergy(self, time_series=False): """ Get the current kinetic energy. Parameters ---------- time_series : bool Whether to return a list of time series records. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The kinetic energy. """ return self.getKineticEnergy(time_series, block=False)
[docs] def getTotalEnergy(self, time_series=False, block="AUTO"): """ Get the total energy. Parameters ---------- time_series : bool Whether to return a list of time series records. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The total energy. """ return self.getRecord( "TOTALENERGY(KJ/MOLE)", time_series, _Units.Energy.kj_per_mol, block )
[docs] def getCurrentTotalEnergy(self, time_series=False): """ Get the current total energy. Parameters ---------- time_series : bool Whether to return a list of time series records. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The total energy. """ return self.getTotalEnergy(time_series, block=False)
[docs] def getTemperature(self, time_series=False, block="AUTO"): """ Get the temperature. Parameters ---------- time_series : bool Whether to return a list of time series records. block : bool Whether to block until the process has finished running. Returns ------- temperature : :class:`Temperature <BioSimSpace.Types.Temperature>` The temperature. """ return self.getRecord( "TEMPERATURE(K)", time_series, _Units.Temperature.kelvin, block )
[docs] def getCurrentTemperature(self, time_series=False): """ Get the current temperature. Parameters ---------- time_series : bool Whether to return a list of time series records. Returns ------- temperature : :class:`Temperature <BioSimSpace.Types.Temperature>` The current temperature. """ return self.getTemperature(time_series, block=False)
[docs] def getVolume(self, time_series=False, block="AUTO"): """ Get the volume. Parameters ---------- time_series : bool Whether to return a list of time series records. block : bool Whether to block until the process has finished running. Returns ------- volume : :class:`Volume <BioSimSpace.Types.Volume>` The volume. """ return self.getRecord( "BOXVOLUME(NM^3)", time_series, _Units.Volume.nanometer3, block )
[docs] def getCurrentVolume(self, time_series=False): """ Get the current volume. Parameters ---------- time_series : bool Whether to return a list of time series records. Returns ------- volume : :class:`Volume <BioSimSpace.Types.Volume>` The volume. """ return self.getVolume(time_series, block=False)
[docs] def stdout(self, n=10): """ Print the last n lines of the stdout buffer. Parameters ---------- n : int The number of lines to print. """ # Note that thermodynamic records, e.g. energy, pressure, temperature, # are redirected to a log file. # Ensure that the number of lines is positive. if n < 0: raise ValueError("The number of lines must be positive!") # Append any new lines to the stdout list. for line in _pygtail.Pygtail(self._stdout_file): self._stdout.append(line.rstrip()) # Get the current number of lines. num_lines = len(self._stdout) # Set the line from which to start printing. if num_lines < n: start = 0 else: start = num_lines - n # Print the lines. for x in range(start, num_lines): print(self._stdout[x])
def _add_config_imports(self): """ Helper function to write the header (import statements) to the OpenMM Python script (config file). """ # We should verify that openmm is available to prevent # difficult-to-debug errors in the run script from BioSimSpace._Utils import _try_import, _assert_imported _openmm = _try_import("openmm") _assert_imported(_openmm) self.addToConfig("from openmm import *") self.addToConfig("from openmm.app import *") self.addToConfig("from openmm.unit import *") self.addToConfig( "from openmm.app.internal.unitcell import reducePeriodicBoxVectors" ) self.addToConfig("import parmed") def _add_config_platform(self): """ Helper function to add platform information to the OpenMM Python script. """ # Set the simulation platform. self.addToConfig("\n# Set the simulation platform.") self.addToConfig(f"platform = Platform.getPlatformByName('{self._platform}')") if self._platform == "CPU": self.addToConfig("properties = {}") elif self._platform == "CUDA": cuda_devices = _os.environ.get("CUDA_VISIBLE_DEVICES") if cuda_devices is None: cuda_devices = "0" _warnings.warn( "'CUDA' platform selected but 'CUDA_VISIBLE_DEVICES' " "environment variable is unset. Defaulting to '0'." ) self.addToConfig("properties = {'CudaDeviceIndex': '%s'}" % cuda_devices) elif self._platform == "OPENCL": opencl_devices = _os.environ.get("OPENCL_VISIBLE_DEVICES") if opencl_devices is None: opencl_devices = "0" _warnings.warn( "'OpenCL' platform selected but 'OPENCL_VISIBLE_DEVICES' " "environment variable is unset. Defaulting to '0'." ) self.addToConfig( "properties = {'OpenCLDeviceIndex': '%s'}" % opencl_devices ) def _add_config_restart(self): """Helper function to check for a restart file and load state information.""" self.addToConfig("\n# Check for a restart file.") self.addToConfig(f"if os.path.isfile('{self._name}.xml'):") self.addToConfig(" is_restart = True") self.addToConfig(f" simulation.loadState('{self._name}.xml')") self.addToConfig(f" if not os.path.isfile('{self._name}.log'):") self.addToConfig(f" raise IOError('Missing log file: {self._name}.log')") self.addToConfig(f" with open('{self._name}.log', 'r') as f:") self.addToConfig(" lines = f.readlines()") self.addToConfig(" last_line = lines[-1].split()") self.addToConfig(" try:") self.addToConfig(" step = int(last_line[0])") self.addToConfig(" except:") self.addToConfig( f" raise IOError('Failed to read current integration step from {self._name}.log')" ) self.addToConfig(" simulation.currentStep = step") self.addToConfig("else:") self.addToConfig(" is_restart = False") if _os.path.isfile(_os.path.join(str(self._work_dir), f"{self._name}.xml")): with open( _os.path.join(str(self._work_dir), f"{self._name}.log"), "r" ) as f: lines = f.readlines() last_line = lines[-1].split() step = int(last_line[0]) return True, step else: return False, 0 def _add_config_monkey_patches(self): """ Helper function to write any monkey-patches to the OpenMM Python script (config file). """ # We monkey-patch the OpenMM DCDFile.writeModel method to avoid writing the # positions of any dummy atoms that are used as restraints to trajectory files. # This avoids the need to delete the dummies from the molecular system on read, # allowing us to make use of System._updateCoordinates which requires that the # number of atoms are consistent between systems. (Deleting the dummies from the # system is slower than not writing them in the first place.) self.addToConfig( "\n# Monkey-patch the DCD.writeModel method to avoid writing dummy-atom positions." ) # Store the original writeModel method. self.addToConfig("writeModel = DCDFile.writeModel") # Create a monkey-patch where we slice the positions list to match the # number of atoms in the topology, then pass this through to the original # writeModel method. self.addToConfig( "def writeModelPatched(self, positions, unitCellDimensions=None, periodicBoxVectors=None):" ) self.addToConfig( " positions = positions[:len(list(self._topology.atoms()))]" ) self.addToConfig(" writeModel(self,") self.addToConfig(" positions,") self.addToConfig(" unitCellDimensions=unitCellDimensions,") self.addToConfig(" periodicBoxVectors=periodicBoxVectors)") # Replace the writeModel method with the monkey-patch. self.addToConfig("DCDFile.writeModel = writeModelPatched") def _add_config_reporters( self, state_interval=100, traj_interval=500, is_restart=False ): """ Helper function to write the reporter (output statements) section to the OpenMM Python script (config file). Parameters ---------- state_interval : int The frequency at which to write state information in integration steps. traj_interval : int The frequency at which to write trajectory frames in integration steps. is_restart : bool Whether the simulation is a restart. """ if not type(state_interval) is int: raise TypeError("'state_interval' must be of type 'int'.") if state_interval <= 0: raise ValueError("'state_interval' must be a positive integer.") if not type(traj_interval) is int: raise TypeError("'traj_interval' must be of type 'int'.") if traj_interval <= 0: raise ValueError("'traj_interval' must be a positive integer.") if not isinstance(is_restart, bool): raise TypeError("'is_restart' must be of type 'bool'.") # Append to a trajectory file every 500 steps. self.addToConfig( f"simulation.reporters.append(DCDReporter('{self._name}.dcd', {traj_interval}, append={is_restart}))" ) # Disable specific state information for minimisation protocols. if isinstance( self._protocol, (_Protocol.Minimisation, _Protocol.ATMMinimisation) ): is_step = False is_time = False is_temperature = False else: is_step = True is_time = True is_temperature = True # Work out the total number of steps. if isinstance( self._protocol, (_Protocol.Minimisation, _Protocol.ATMMinimisation) ): total_steps = 1 else: total_steps = _math.ceil( self._protocol.getRunTime() / self._protocol.getTimeStep() ) # Write state information to file every 100 steps. self.addToConfig(f"log_file = open('{self._name}.log', 'a')") self.addToConfig(f"simulation.reporters.append(StateDataReporter(log_file,") self.addToConfig( f" {state_interval}," ) self.addToConfig( f" step={is_step}," ) self.addToConfig( f" time={is_time}," ) self.addToConfig( " potentialEnergy=True," ) self.addToConfig( " kineticEnergy=True," ) self.addToConfig( " totalEnergy=True," ) self.addToConfig( f" temperature={is_temperature}," ) self.addToConfig(" volume=True,") self.addToConfig( f" totalSteps={total_steps}," ) self.addToConfig(" speed=True,") self.addToConfig( " remainingTime=True," ) self.addToConfig( " separator=' '))" ) def _add_config_restraints(self): """Helper function to add position restraints to the OpenMM script (config file).""" # Add position restraints. This uses the approach from: # https://github.com/openmm/openmm/issues/2262#issuecomment-464157489 # Here zero-mass dummy atoms are bonded to the restrained atoms to avoid # issues with position rescaling during barostat updates. restraint = self._protocol.getRestraint() if restraint is not None: # Search for the atoms to restrain by keyword. if isinstance(restraint, str): restrained_atoms = self._reference_system.getRestraintAtoms(restraint) # Use the user-defined list of indices. else: restrained_atoms = restraint self.addToConfig( f"ref_prm = parmed.load_file('{self._name}.prm7', '{self._name}_ref.rst7')" ) # Get the force constant in units of kJ_per_mol/nanometer**2 force_constant = self._protocol.getForceConstant()._sire_unit force_constant = force_constant.to( _SireUnits.kJ_per_mol / _SireUnits.nanometer2 ) self.addToConfig( "\n# Restrain the position of atoms using zero-mass dummy atoms." ) self.addToConfig("restraint = HarmonicBondForce()") self.addToConfig("restraint.setUsesPeriodicBoundaryConditions(True)") self.addToConfig("system.addForce(restraint)") self.addToConfig( "nonbonded = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]" ) self.addToConfig("dummy_indices = []") self.addToConfig("positions = ref_prm.positions") self.addToConfig(f"restrained_atoms = {restrained_atoms}") self.addToConfig("for i in restrained_atoms:") self.addToConfig(" j = system.addParticle(0)") self.addToConfig(" nonbonded.addParticle(0, 1, 0)") self.addToConfig(" nonbonded.addException(i, j, 0, 1, 0)") self.addToConfig( f" restraint.addBond(i, j, 0*nanometers, {force_constant}*kilojoules_per_mole/nanometer**2)" ) self.addToConfig(" dummy_indices.append(j)") self.addToConfig(" positions.append(positions[i])") def _update_stdout_dict(self): """Update the dictionary of thermodynamic records.""" # Exit if log file hasn't been created. if not _os.path.isfile(self._log_file): return # A list of the new record lines. lines = [] # Append any new lines. for line in _pygtail.Pygtail(self._log_file): lines.append(line) # Append any new records to the stdout dictionary. for line in lines: # Strip leading/trailing whitespace. line = line.strip() # This is the header record. if line[0] == "#": if self._header is None: # Store the header. self._header = line else: # If this is a restart, make sure the information in the file # is consistent. if line != self._header: raise _IncompatibleError( "Mismatch in the log file header.! " f"Original header: '{self._header}', " f"Current header: '{line}'" ) # Work out what records are in the file and the separator # that is used. While we use a standard format, this makes # sure that we can still parse the log file if the user # happens to have changed the formatting. # A tally counter for the number of quotes that we've seen # in the line so far. num_quotes = 0 # Initialise the separator. sep = "" # Loop over the characters in the line. for c in line: # Increment the number of quotes. if c == '"': num_quotes += 1 # This is the second quote we've seen, start adding # characters to the separator. if num_quotes == 2: sep += c # Break when we've reached the next quote. elif num_quotes == 3: break # The separator includes a leading " character, so delete it. self._record_separator = sep[1:] # Now split the line on the separator to work out the records. # We ignore the first character since it is a comment. # Here we use the full separator, i.e. including the both quotes, # so that we can correctly split record names with spaces in them. records = line[1:].split(sep + '"') # Store the number of records. num_records = len(records) # Map each record string to its position in the array (column order). for idx, record in enumerate(records): # Strip the extra quotes from the record. if idx == 0: record = record[1:] elif idx == num_records - 1: record = record[:-1] # Map the index to the record. Store records in upper # case without whitespace to help catch typos from # the user. self._record_mapping[idx] = record.upper().replace(" ", "") # Extract the records and add them to the dictionary. else: # Split the line on the separator. records = line.split(self._record_separator) # Add each record to the appropriate key in the MultiDict. for idx, record in enumerate(records): # Get the key for this record. key = self._record_mapping[idx] # Update the record dictionary for this key. self._stdout_dict[key] = record def _get_stdout_record(self, key, time_series=False, unit=None): """ Helper function to get a stdout record from the dictionary. Parameters ---------- key : str The record key. time_series : bool Whether to return a time series of records. unit : BioSimSpace.Types._type.Type The unit to convert the record to. Returns ------- record : The matching stdout record. """ # No data! if len(self._stdout_dict) == 0: return None if not isinstance(time_series, bool): _warnings.warn("Non-boolean time-series flag. Defaulting to False!") time_series = False # Validate the unit. if unit is not None: if not isinstance(unit, _Type): raise TypeError("'unit' must be of type 'BioSimSpace.Types'") # Return the list of dictionary values. if time_series: try: if unit is None: return [float(x) for x in self._stdout_dict[key]] else: return [ (float(x) * unit)._to_default_unit() for x in self._stdout_dict[key] ] except KeyError: return None # Return the most recent dictionary value. else: try: if unit is None: return float(self._stdout_dict[key][-1]) else: return (float(self._stdout_dict[key][-1]) * unit)._to_default_unit() except KeyError: return None