Source code for BioSimSpace.Process._amber

######################################################################
# 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.
#
# BioSimSpace 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 using AMBER."""

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

__all__ = ["Amber"]

from .._Utils import _try_import

_pygtail = _try_import("pygtail")

import os as _os
import re as _re
import time as _time
import shutil as _shutil
import tempfile as _tempfile
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 Mol as _SireMol

from .. import _amber_home, _isVerbose
from ..Align._squash import _squash, _unsquash
from .._Config import Amber as _AmberConfig
from .._Exceptions import IncompatibleError as _IncompatibleError
from .._Exceptions import MissingSoftwareError as _MissingSoftwareError
from ..Protocol._free_energy_mixin import _FreeEnergyMixin
from ..Protocol._position_restraint_mixin import _PositionRestraintMixin
from .._SireWrappers import System as _System
from ..Types._type import Type as _Type

from .. import IO as _IO
from .. import Protocol as _Protocol
from .. import Trajectory as _Trajectory
from .. import Units as _Units
from .. import _Utils

from . import _process

from ._plumed import Plumed as _Plumed


[docs] class Amber(_process.Process): """A class for running simulations using AMBER."""
[docs] def __init__( self, system, protocol, reference_system=None, explicit_dummies=False, exe=None, is_gpu=False, name="amber", work_dir=None, seed=None, extra_options={}, extra_lines=[], extra_args={}, property_map={}, **kwargs, ): """ Constructor. Parameters ---------- system : :class:`System <BioSimSpace._SireWrappers.System>` The molecular system. protocol : :class:`Protocol <BioSimSpace.Protocol>` The protocol for the AMBER 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. explicit_dummies : bool Whether to keep dummy atoms explicit at alchemical end states, or remove them. exe : str The full path to the AMBER executable. is_gpu : bool Whether to use the GPU accelerated version of AMBER. name : str The name of the process. work_dir : The working directory for the process. seed : int A random number seed. extra_options : dict A dictionary containing extra options. Overrides the defaults generated by the protocol. extra_lines : [str] A list of extra lines to put at the end of the configuration file. extra_args : dict A dictionary of extra command-line arguments to pass to the AMBER executable. 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, extra_options=extra_options, extra_lines=extra_lines, extra_args=extra_args, property_map=property_map, ) # Set the package name. self._package_name = "AMBER" # This process can generate trajectory data. self._has_trajectory = True if not isinstance(is_gpu, bool): raise TypeError("'is_gpu' must be of type 'bool'") # Check whether this is a vacuum simulation. self._is_vacuum = not ( _AmberConfig.hasBox(self._system, self._property_map) or _AmberConfig.hasWater(self._system) ) # Flag to indicate whether the original system has a box. self._has_box = _AmberConfig.hasBox(self._system, self._property_map) # If the path to the executable wasn't specified, then search # for it in AMBERHOME and the PATH. if exe is None: if isinstance(protocol, _FreeEnergyMixin): is_free_energy = True else: is_free_energy = False self._exe = _find_exe( is_gpu=is_gpu, is_free_energy=is_free_energy, is_vacuum=self._is_vacuum ) else: # Make sure executable exists. if _os.path.isfile(exe): self._exe = exe else: raise IOError("AMBER executable doesn't exist: '%s'" % exe) # Is this a CUDA enabled version of AMBER? if "cuda" in self._exe.lower(): self._is_pmemd_cuda = True self._is_pmemd = False else: self._is_pmemd_cuda = False if "pmemd" in self._exe.lower(): self._is_pmemd = True else: self._is_pmemd = False if not isinstance(explicit_dummies, bool): raise TypeError("'explicit_dummies' must be of type 'bool'") self._explicit_dummies = explicit_dummies # Initialise the energy dictionary and header. self._stdout_dict = _process._MultiDict() # Initialise dictionaries to hold stdout records for all possible # regions. For regular simulations there will be one, for free-energy # simulations there can be up to four, i.e. one for each of the TI regions # and one for the soft-core part of the system in each region, if present. # The order of the dictionaries is: # - TI region 1 # - TI region 1 (soft-core part) # - TI region 2 # - TI region 2 (soft-core part) self._stdout_dict = [ _process._MultiDict(), _process._MultiDict(), _process._MultiDict(), _process._MultiDict(), ] # Initialise mappings between "universal" stdout keys, and the actual # record key used for the different regions (and soft-core parts) from # in the AMBER output. Ordering is the same as for the stdout_dicts above. self._stdout_key = [{}, {}, {}, {}] # Flag for the current record region in the AMBER output file. self._current_region = 0 # Initialise log file parsing flags. self._has_results = False self._finished_results = False self._is_header = False # The names of the input 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}.nc") # Set the path for the AMBER configuration file. self._config_file = _os.path.join(str(self._work_dir), f"{name}.cfg") # 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) # Now set up the working directory for the process. self._setup(**kwargs)
def _setup(self, **kwargs): """Setup the input files and working directory ready for simulation.""" # Create the input files... # 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 ) # Create the squashed system. if isinstance(self._protocol, _FreeEnergyMixin): # Check that the system contains a perturbable molecule. if self._system.nPerturbableMolecules() == 0: raise ValueError( "'BioSimSpace.Protocol.FreeEnergy' requires a " "perturbable molecule!" ) # 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: # Get the existing box information. box, _ = system.getBox() # We need to add a box. if box is None: from ..Box import cubic as _cubic from ..Units.Length import angstrom as _angstrom # Get the bounding box of the system. box_min, box_max = system.getAxisAlignedBoundingBox() # Work out the box size from the difference in the coordinates. box_size = [y - x for x, y in zip(box_min, box_max)] # Work out the size of the box assuming an 8 Angstrom non-bonded cutoff. padding = 8 * _angstrom box_length = max(box_size) + padding # Small box fix. Should be patched in future versions of pmemd.cuda. if box_length < 30 * _angstrom: box_length = 30 * _angstrom # Set the simulation box. system.setBox(*_cubic(box_length)) # Apply SOMD1 compatibility to the perturbation. if ( "somd1_compatibility" in kwargs and kwargs.get("somd1_compatibility") is True ): from ._somd import _somd1_compatibility system = _somd1_compatibility(system) system, self._mapping = _squash( system, explicit_dummies=self._explicit_dummies ) self._squashed_system = system else: # Check for perturbable molecules and convert to the chosen end state. system = self._checkPerturbable(system) # 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 file for position restraints. try: file = _os.path.splitext(self._ref_file)[0] _IO.saveMolecules( file, self._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 # Generate the AMBER configuration file. # 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 AMBER configuration file strings.""" extra_options = self._extra_options.copy() extra_lines = self._extra_lines.copy() # Set the random number seed. if self._seed is None: extra_options["ig"] = -1 else: extra_options["ig"] = self._seed # Add configuration variables for a metadynamics simulation. if isinstance(self._protocol, (_Protocol.Metadynamics, _Protocol.Steering)): extra_options["plumed"] = 1 extra_options["plumedfile"] = "'plumed.dat'" # Create the PLUMED input file and copy auxiliary files to the working directory. self._plumed = _Plumed(str(self._work_dir)) plumed_config, auxiliary_files = self._plumed.createConfig( self._system, self._protocol, self._property_map ) self._setPlumedConfig(plumed_config) if auxiliary_files is not None: for file in auxiliary_files: file_name = _os.path.basename(file) _shutil.copyfile( file, _os.path.join(str(self._work_dir), file_name) ) self._input_files.append(self._plumed_config_file) # Expose the PLUMED specific member functions. setattr(self, "getPlumedConfig", self._getPlumedConfig) setattr(self, "getPlumedConfigFile", self._getPlumedConfigFile) setattr(self, "setPlumedConfig", self._setPlumedConfig) setattr(self, "getFreeEnergy", self._getFreeEnergy) setattr(self, "getCollectiveVariable", self._getCollectiveVariable) setattr(self, "sampleConfigurations", self._sampleConfigurations) setattr(self, "getTime", self._getTime) # Instantiate the AMBER configuration generator. amber_config = _AmberConfig(self._system, self._protocol) # Create the configuration. self.setConfig( amber_config.createConfig( is_pmemd=self._is_pmemd, is_pmemd_cuda=self._is_pmemd_cuda, explicit_dummies=self._explicit_dummies, extra_options=extra_options, extra_lines=extra_lines, ) ) # Flag that this isn't a custom protocol. if not self._extra_options and not self._extra_lines: self._protocol._setCustomised(False) def _generate_args(self): """Generate the dictionary of command-line arguments.""" # Clear the existing arguments. self.clearArgs() # Add the default arguments. self.setArg("-O", True) # Overwrite. self.setArg("-i", "%s.cfg" % self._name) # Input file. self.setArg("-p", "%s.prm7" % self._name) # Topology file. self.setArg("-c", "%s.rst7" % self._name) # Coordinate file. self.setArg("-o", "%s.out" % self._name) # Redirect stdout to file. self.setArg("-r", "%s.crd" % self._name) # Restart file. self.setArg("-inf", "%s.nrg" % self._name) # Energy info file. # Skip if the user has passed a custom protocol. if not isinstance(self._protocol, _Protocol.Custom): # Append a reference file if a position restraint is specified. if isinstance(self._protocol, _PositionRestraintMixin): if self._protocol.getRestraint() is not None: self.setArg("-ref", "%s_ref.rst7" % self._name) # Append a trajectory file if this anything other than a minimisation. if not isinstance(self._protocol, _Protocol.Minimisation): self.setArg("-x", "%s.nc" % self._name) # Add the extra arguments. for key, value in self._extra_args.items(): self.setArg(key, value)
[docs] def start(self): """ Start the AMBER process. Returns ------- process : :class:`Process.Amber <BioSimSpace.Process.Amber>` The process object. """ # The process is currently queued. if self.isQueued(): return # Process is already running. if self._process is not None: if self._process.isRunning(): return # Run the process in the working directory. with _Utils.cd(self._work_dir): # Create the arguments string list. args = self.getArgStringList() # Write the command-line process to a README.txt file. with open("README.txt", "w") as file: # Set the command-line string. self._command = "%s " % self._exe + self.getArgString() # Write the command to file. file.write("# AMBER was run with the following command:\n") file.write("%s\n" % self._command) # Start the timer. self._timer = _timeit.default_timer() # Start the simulation. Pass a null string for the stdout file # since we've explicitly redirected AMBER output to file since # pmemd doesn't write to standard output. self._process = _SireBase.Process.run( self._exe, args, "", "%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() # Warn the user if the process has exited with an error. if self.isError(): _warnings.warn("The process exited with an error!") # Create the name of the restart CRD file. restart = _os.path.join(str(self._work_dir), "%s.crd" % self._name) # Check that the file exists. if _os.path.isfile(restart): # 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 # Copy the restart file to a temporary location. Sire streams from # binary files with the same path, so we need to ensure that a new # stream is created each time. with _tempfile.TemporaryDirectory() as tmp_dir: tmp_file = f"{tmp_dir}/{self._name}.crd" _shutil.copyfile(restart, tmp_file) # Create a new molecular system from the restart file. new_system = _System( _SireIO.MoleculeParser.read( [tmp_file, self._top_file], self._property_map ) ) # Create a copy of the existing system object. old_system = self._system.copy() if isinstance(self._protocol, _FreeEnergyMixin): # Udpate the coordinates and velocities and return a mapping between # the molecule indices in the two systems. mapping = { _SireMol.MolIdx(x): _SireMol.MolIdx(x) for x in range(0, self._squashed_system.nMolecules()) } ( self._squashed_system._sire_object, _, ) = _SireIO.updateCoordinatesAndVelocities( self._squashed_system._sire_object, new_system._sire_object, mapping, is_lambda1, self._property_map, self._property_map, ) # Update the unsquashed system based on the updated squashed system. old_system = _unsquash( old_system, self._squashed_system, self._mapping, explicit_dummies=self._explicit_dummies, ) else: # 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 self._has_box: 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: 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:`Trajectory <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!") try: return _Trajectory.Trajectory(process=self, backend=backend) except: return None
[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() ) if index < 0 or index > max_index: raise ValueError(f"'index' must be in range [0, {max_index}].") 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 # Get the latest trajectory frame. new_system = _Trajectory.getFrame(self._traj_file, self._top_file, index) # Create a copy of the existing system object. 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 self._has_box: 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 getRecordKey(self, record, region=0, soft_core=False): """ Parameters ---------- record : str The record used in the AMBER standard output, e.g. 'TEMP(K)'. Please consult the current AMBER manual for details: https://ambermd.org/Manuals.php region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- key : str The universal record key that can be used with getRecord. """ # Validate the record string. if not isinstance(record, str): raise TypeError("'record' must be of type 'str'") # Validate the region. if not isinstance(region, int): raise TypeError("'region' must be of type 'int'") else: if region < 0 or region > 1: raise ValueError("'region' must be in range [0, 1]") # Validate the soft-core flag. if not isinstance(soft_core, bool): raise TypeError("'soft_core' must be of type 'bool'.") # Convert to the full index. idx = 2 * region + int(soft_core) # Strip whitespace from the beginning and end of the record and convert # to upper case. cleaned_record = record.strip().upper() # Make sure the record exists in the key mapping. if not cleaned_record in self._stdout_key[idx].values(): raise ValueError(f"No key found for record '{record}'") return list(self._stdout_key[idx].keys())[ list(self._stdout_key[idx].values()).index(cleaned_record) ]
[docs] def getRecord( self, key, time_series=False, unit=None, region=0, soft_core=False, block="AUTO" ): """ Get a record from the stdout dictionary. Parameters ---------- key : str A universal record key based on the key used in the AMBER standard output. Use 'getRecordKey(record)` to generate the key. The records are those used in the AMBER standard output, e.g. 'TEMP(K)'. Please consult the current AMBER manual for details: https://ambermd.org/Manuals.php time_series : bool Whether to return a list of time series records. unit : :class:`Unit <BioSimSpace.Units>` The unit to convert the record to. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. 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!") return self._get_stdout_record( key.strip().upper(), time_series=time_series, unit=unit, region=region, soft_core=soft_core, )
[docs] def getCurrentRecord( self, key, time_series=False, unit=None, region=0, soft_core=False ): """ Get a current record from the stdout dictionary. Parameters ---------- key : str A universal record key based on the key used in the AMBER standard output. Use 'getRecordKey(record)` to generate the key. The records are those used in the AMBER standard output, e.g. 'TEMP(K)'. Please consult the current AMBER manual for details: https://ambermd.org/Manuals.php time_series : bool Whether to return a list of time series records. unit : :class:`Unit <BioSimSpace.Units>` The unit to convert the record to. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. 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!") return self._get_stdout_record( key.strip().upper(), time_series=time_series, unit=unit, region=region, soft_core=soft_core, )
[docs] def getRecords(self, region=0, soft_core=False, block="AUTO"): """ Return the dictionary of stdout time-series records. Parameters ---------- region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. 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. """ # Validate the region. if not isinstance(region, int): raise TypeError("'region' must be of type 'int'") else: if region < 0 or region > 1: raise ValueError("'region' must be in range [0, 1]") # Validate the soft-core flag. if not isinstance(soft_core, bool): raise TypeError("'soft_core' must be of type 'bool'.") # Convert to the full index, region + soft_core. idx = 2 * region + int(soft_core) # 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.stdout(0) return self._stdout_dict[idx].copy()
[docs] def getCurrentRecords(self, region=0, soft_core=False): """ Return the current dictionary of stdout time-series records. Parameters ---------- region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- records : :class:`MultiDict <BioSimSpace.Process._process._MultiDict>` The dictionary of time-series records. """ return self.getRecords(region=region, soft_core=soft_core, block=False)
[docs] def getTime(self, time_series=False, region=0, soft_core=False, block="AUTO"): """ Get the simulation time. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- time : :class:`Time <BioSimSpace.Types.Time>` The current simulation time in nanoseconds. """ # No time records for minimisation protocols. if isinstance(self._protocol, _Protocol.Minimisation): return None # Get the list of time steps. time_steps = self.getRecord( "TIME(PS)", time_series=time_series, unit=None, region=region, soft_core=soft_core, block=block, ) # Convert from picoseconds to nanoseconds. if time_steps is not None: if time_series: return [ (x * _Units.Time.picosecond)._to_default_unit() for x in time_steps ] else: return (time_steps * _Units.Time.picosecond)._to_default_unit()
[docs] def getCurrentTime(self, time_series=False, region=0, soft_core=False): """ Get the current simulation time. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- time : :class:`Time <BioSimSpace.Types.Time>` The current simulation time in nanoseconds. """ return self.getTime( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getStep(self, time_series=False, region=0, soft_core=False, block="AUTO"): """ Get the number of integration steps. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- step : int The current number of integration steps. """ return self.getRecord( "NSTEP", time_series=time_series, unit=None, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentStep(self, time_series=False, region=0, soft_core=False): """ Get the current number of integration steps. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- step : int The current number of integration steps. """ return self.getStep( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getBondEnergy(self, time_series=False, region=0, soft_core=False, block="AUTO"): """ Get the bond energy. Parameters ---------- time_series : bool Whether to return a list of time series records. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The bond energy. """ return self.getRecord( "BOND", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentBondEnergy(self, time_series=False, region=0, soft_core=False): """ Get the current bond energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The bond energy. """ return self.getBondEnergy( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getAngleEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the angle energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The angle energy. """ return self.getRecord( "ANGLE", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentAngleEnergy(self, time_series=False, region=0, soft_core=False): """ Get the current angle energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The angle energy. """ return self.getAngleEnergy( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getDihedralEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the total dihedral energy (proper + improper). Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The total dihedral energy. """ return self.getRecord( "DIHED", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentDihedralEnergy(self, time_series=False, region=0, soft_core=False): """ Get the current total dihedral energy (proper + improper). Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The total dihedral energy. """ return self.getDihedralEnergy( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getElectrostaticEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the electrostatic energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The electrostatic energy. """ return self.getRecord( "EEL", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentElectrostaticEnergy( self, time_series=False, region=0, soft_core=False ): """ Get the current dihedral energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The electrostatic energy. """ return self.getElectrostaticEnergy( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getElectrostaticEnergy14( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the electrostatic energy between atoms 1 and 4. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The electrostatic energy between atoms 1 and 4. """ return self.getRecord( "14EEL", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentElectrostaticEnergy14( self, time_series=False, region=0, soft_core=False ): """ Get the current electrostatic energy between atoms 1 and 4. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The electrostatic energy between atoms 1 and 4. """ return self.getElectrostaticEnergy14( time_series=time_series, region=region, soft_core=False, block=False )
[docs] def getVanDerWaalsEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the Van der Vaals energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The Van der Vaals energy. """ return self.getRecord( "VDW", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentVanDerWaalsEnergy(self, time_series=False, region=0, soft_core=False): """ Get the current Van der Vaals energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The Van der Vaals energy. """ return self.getVanDerWaalsEnergy( time_series=time_series, block=False, region=region, soft_core=soft_core )
[docs] def getVanDerWaalsEnergy14( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the Van der Vaals energy between atoms 1 and 4. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The Van der Vaals energy between atoms 1 and 4. """ return self.getRecord( "14VDW", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentVanDerWaalsEnergy14( self, time_series=False, region=0, soft_core=False ): """ Get the current Van der Vaals energy between atoms 1 and 4. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The Van der Vaals energy between atoms 1 and 4. """ return self.getVanDerWaalsEnergy( time_series=time_series, block=False, region=region, soft_core=soft_core )
[docs] def getHydrogenBondEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the hydrogen bond energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The hydrogen bond energy. """ return self.getRecord( "EHBOND", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentHydrogenBondEnergy( self, time_series=False, region=0, soft_core=False ): """ Get the current hydrogen bond energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The hydrogen bond energy. """ return self.getHydrogenBondEnergy( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getRestraintEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the restraint energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The restraint energy. """ return self.getRecord( "RESTRAINT", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentRestraintEnergy(self, time_series=False, region=0, soft_core=False): """ Get the current restraint energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The restraint energy. """ return self.getRestraintEnergy( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getPotentialEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the potential energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The potential energy. """ return self.getRecord( "EPTOT", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentPotentialEnergy(self, time_series=False, region=0, soft_core=False): """ Get the current potential energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The potential energy. """ return self.getPotentialEnergy( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getKineticEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the kinetic energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The kinetic energy. """ return self.getRecord( "EKTOT", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentKineticEnergy(self, time_series=False, region=0, soft_core=False): """ Get the current kinetic energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The kinetic energy. """ return self.getKineticEnergy( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getNonBondedEnergy14( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the non-bonded energy between atoms 1 and 4. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The non-bonded energy between atoms 1 and 4. """ return self.getRecord( "14NB", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentNonBondedEnergy14(self, time_series=False, region=0, soft_core=False): """ Get the current non-bonded energy between atoms 1 and 4. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The non-bonded energy between atoms 1 and 4. """ return self.getNonBondedEnergy14( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getTotalEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the total energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The total energy. """ if not isinstance(region, int): raise TypeError("'region' must be of type 'int'") else: if region < 0 or region > 1: raise ValueError("'region' must be in range [0, 1]") # Validate the soft-core flag. if not isinstance(soft_core, bool): raise TypeError("'soft_core' must be of type 'bool'.") # Convert to the full index, region + soft_core. idx = 2 * region + int(soft_core) if isinstance(self._protocol, _Protocol.Minimisation) and not soft_core: return self.getRecord( "ENERGY", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, ) else: return self.getRecord( "ETOT", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentTotalEnergy(self, time_series=False, region=0, soft_core=False): """ Get the current total energy. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The total energy. """ return self.getTotalEnergy( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getCentreOfMassKineticEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the kinetic energy of the centre of mass in translation. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The centre of mass kinetic energy. """ return self.getRecord( "EKCMT", time_series=time_series, unit=_Units.Energy.kcal_per_mol, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentCentreOfMassKineticEnergy( self, time_series=False, region=0, soft_core=False ): """ Get the current kinetic energy of the centre of mass in translation. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- energy : :class:`Energy <BioSimSpace.Types.Energy>` The centre of mass kinetic energy. """ return self.getCentreOfMassKineticEnergy( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getVirial(self, time_series=False, region=0, soft_core=False, block="AUTO"): """ Get the virial. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- virial : float The virial. """ return self.getRecord( "VIRIAL", time_series=time_series, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentVirial(self, time_series=False, region=0, soft_core=False): """ Get the current virial. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- virial : float The virial. """ return self.getVirial( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getTemperature( self, time_series=False, region=0, soft_core=False, block="AUTO" ): """ Get the temperature. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- temperature : :class:`Temperature <BioSimSpace.Types.Temperature>` The temperature. """ return self.getRecord( "TEMP(K)", time_series=time_series, unit=_Units.Temperature.kelvin, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentTemperature(self, time_series=False, region=0, soft_core=False): """ Get the current temperature. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- temperature : :class:`Temperature <BioSimSpace.Types.Temperature>` The temperature. """ return self.getTemperature( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getPressure(self, time_series=False, region=0, soft_core=False, block="AUTO"): """ Get the pressure. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- pressure : :class:`Pressure <BioSimSpace.Types.Pressure>` The pressure. """ return self.getRecord( "PRESS", time_series=time_series, unit=_Units.Pressure.bar, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentPressure(self, time_series=False, region=0, soft_core=False): """ Get the current pressure. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- pressure : :class:`Pressure <BioSimSpace.Types.Pressure>` The pressure. """ return self.getPressure( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getVolume(self, time_series=False, region=0, soft_core=False, block="AUTO"): """ Get the volume. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- volume : :class:`Volume <BioSimSpace.Types.Volume>` The volume. """ return self.getRecord( "VOLUME", time_series=time_series, unit=_Units.Volume.angstrom3, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentVolume(self, time_series=False, region=0, soft_core=False): """ Get the current volume. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- volume : :class:`Volume <BioSimSpace.Types.Volume>` The volume. """ return self.getVolume( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getDensity(self, time_series=False, region=0, soft_core=False, block="AUTO"): """ Get the density. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- density : float The density. """ return self.getRecord( "DENSITY", time_series=time_series, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentDensity(self, time_series=False, region=0, soft_core=False): """ Get the current density. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- density : float The density. """ return self.getDensity( time_series=time_series, region=region, soft_core=soft_core, block=False )
[docs] def getDVDL(self, time_series=False, region=0, soft_core=False, block="AUTO"): """ Get the gradient of the total energy with respect to lambda. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. block : bool Whether to block until the process has finished running. Returns ------- dv_dl : float The gradient of the total energy with respect to lambda. """ return self.getRecord( "DVDL", time_series=time_series, region=region, soft_core=soft_core, block=block, )
[docs] def getCurrentDVDL(self, time_series=False, region=0, soft_core=False): """ Get the current gradient of the total energy with respect to lambda. Parameters ---------- time_series : bool Whether to return a list of time series records. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- dv_dl : float The current gradient of the total energy with respect to lambda. """ return self.getDVDL( time_series=time_series, region=region, soft_core=soft_core, 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. """ # Ensure that the number of lines is positive. if n < 0: raise ValueError("The number of lines must be positive!") # Flag that this isn't a header line. self._is_header = False # Append any new lines to the stdout list. for line in _pygtail.Pygtail(self._stdout_file): self._stdout.append(line.rstrip()) line = line.strip() # Swap dictionary based on the protocol and the degre of freedom to # which the next block of records correspond. if isinstance(self._protocol, _FreeEnergyMixin): if "TI region 1" in line: self._current_region = 0 elif "TI region 2" in line: self._current_region = 2 elif "Softcore part" in line and self._current_region == 0: self._current_region = 1 elif "Softcore part" in line and self._current_region == 2: self._current_region = 3 elif "Detailed TI info" in line: # This flags that we should skip records until the start of # the next set for the first TI region. self._current_region = 4 # Default stdout dictionary. else: self._current_region = 0 # Continue if we are ignoring this record block. if self._current_region == 4: continue stdout_dict = self._stdout_dict[self._current_region] stdout_key = self._stdout_key[self._current_region] # Skip empty lines and summary reports. if len(line) > 0 and line[0] != "|" and line[0] != "-": # Skip EAMBER records. if "EAMBER (non-restraint)" in line: continue # Flag that we've started recording results. elif not self._has_results and line.startswith("NSTEP"): self._has_results = True self._finished_results = False # Flag that we've finished recording results. elif "A V E R A G E S" in line: self._finished_results = True # Parse the results. if self._has_results and not self._finished_results: # The first line of output has different formatting for minimisation protocols. if isinstance(self._protocol, _Protocol.Minimisation): # No equals sign in the line. if "NSTEP" in line and "=" not in line: # Split the line using whitespace. data = line.upper().split() # If we find a header, jump to the top of the loop. if len(data) > 0: if data[0] == "NSTEP": self._is_header = True continue # Process the header record. if self._is_header: # Split the line using whitespace. data = line.upper().split() # The file hasn't been updated. if ( "NSTEP" in stdout_dict and data[0] == stdout_dict["NSTEP"][-1] ): self._finished_results = True continue # Add the timestep and energy records to the dictionary. stdout_dict["NSTEP"] = data[0] stdout_dict["ENERGY"] = data[1] # Add the keys to the mapping stdout_key["NSTEP"] = "NSTEP" stdout_key["ENERGY"] = "ENERGY" # Turn off the header flag now that the data has been recorded. self._is_header = False # All other records are formatted as RECORD = VALUE. # Use a regex search to split the line into record names and values. records = _re.findall( r"([SC_]*[EEL_]*[RES_]*[VDW_]*\d*\-*\d*\s*[A-Z/]+\(*[A-Z]*\)*)\s*=\s*(\-*\d+\.?\d*|\**)", line.upper(), ) # Append each record to the dictionary. for key, value in records: # Strip whitespace from beginning and end. key = key.strip() # Format key so it can be re-used for records corresponding to # different regions, which use different abbreviations. universal_key = ( key.replace("SC_", "") .replace(" ", "") .replace("-", "") .replace("EELEC", "EEL") .replace("VDWAALS", "VDW") ) # Handle missing values, which will appear as asterisks, e.g. # PRESS=******** try: tmp = float(value) except: value = None # Store the record using the original key. stdout_dict[key] = value # Map the universal key to the original. stdout_key[universal_key] = key # 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])
[docs] def kill(self): """Kill the running process.""" # Kill the process. if not self._process is None and self._process.isRunning(): self._process.kill()
def _get_stdout_record( self, key, time_series=False, unit=None, region=0, soft_core=False ): """ Helper function to get a stdout record from the dictionary. Parameters ---------- key : str The universal record key. time_series : bool Whether to return a time series of records. unit : :class:`Type <BioSimSpace.Types._type.Type>` The unit to convert the record to. region : int The region to which the record corresponds. There will only be more than one region for FreeEnergy protocols, where 1 indicates the second TI region. soft_core : bool Whether to get the record for the soft-core part of the system for the chosen region. Returns ------- record : The matching stdout record. """ # Update the standard output dictionary. self.stdout(0) # 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'") # Validate the region. if not isinstance(region, int): raise TypeError("'region' must be of type 'int'") else: if region < 0 or region > 1: raise ValueError("'region' must be in range [0, 1]") # Validate the soft-core flag. if not isinstance(soft_core, bool): raise TypeError("'soft_core' must be of type 'bool'.") # Convert to the full index, region + soft_core. idx = 2 * region + int(soft_core) # Extract the dictionary of stdout records for the specified region and soft-core flag. stdout_dict = self._stdout_dict[idx] # Map the universal key to the original key used for this region. try: key = self._stdout_key[idx][key] except: return None # Return the list of dictionary values. if time_series: try: if key == "NSTEP": return [int(x) for x in stdout_dict[key]] else: if unit is None: return [float(x) if x else None for x in stdout_dict[key]] else: return [ (float(x) * unit)._to_default_unit() if x else None for x in stdout_dict[key] ] except KeyError: return None # Return the most recent dictionary value. else: try: if key == "NSTEP": return int(stdout_dict[key][-1]) else: if unit is None: try: return float(stdout_dict[key][-1]) except: return None else: try: return ( float(stdout_dict[key][-1]) * unit )._to_default_unit() except: return None except KeyError: return None
def _find_exe(is_gpu=False, is_free_energy=False, is_vacuum=False): """ Helper function to search for an AMBER executable. Parameters ---------- is_gpu : bool Whether to search for a GPU-enabled executable. is_free_energy : bool Whether this is a free energy simulation. is_vacuum : bool Whether this is a vacuum simulation. Returns ------- exe : str The path to the executable. """ if not isinstance(is_gpu, bool): raise TypeError("'is_gpu' must be of type 'bool'.") if not isinstance(is_free_energy, bool): raise TypeError("'is_free_energy' must be of type 'bool'.") if not isinstance(is_vacuum, bool): raise TypeError("'is_vacuum' must be of type 'bool'.") if is_gpu: targets = ["pmemd.cuda"] else: if is_free_energy: targets = ["pmemd"] else: targets = ["pmemd", "sander"] # Search for the executable. import os as _os import pathlib as _pathlib from glob import glob as _glob # Get the current path. path = _os.environ["PATH"].split(_os.pathsep) # If AMBERHOME is set, then prepend to the path. if _amber_home is not None: path = [_amber_home + "/bin"] + path # Helper function to check whether a file is executable. def is_exe(fpath): return _os.path.isfile(fpath) and _os.access(fpath, _os.X_OK) # Loop over each directory in the path and search for the executable. for p in path: # Loop over each target. for t in targets: # Glob for the executable. results = _glob(f"{t}*", root_dir=p) # If we find a match, check that it's executable and return the path. # Note that this returns the first match, not the best match. If a # user requires a specific version of the executable, they should # order their path accordingly, or use the exe keyword argument. if results: for exe in results: exe = _pathlib.Path(p) / exe if is_exe(exe): return str(exe) msg = ( "'BioSimSpace.Process.Amber' is not supported. " "Unable to find AMBER executable in AMBERHOME or PATH. " "Please install AMBER (http://ambermd.org)." ) if is_free_energy: msg += " Free energy simulations require 'pmemd' or 'pmemd.cuda'." # If we don't find the executable, raise an error. raise _MissingSoftwareError(msg)