######################################################################
# BioSimSpace: Making biomolecular simulation a breeze!
#
# Copyright: 2017-2025
#
# 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.
This option is provided for legacy support of alchemical free energy calculations
using the old PMEMD CPU implementation. The default is False, which should be
used for any recent AMBER version or for GPU accelerated PMEMD.
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()
reference_system = self._reference_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)
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!"
)
# Make sure the protocol is valid.
if self._protocol.getPerturbationType() != "full":
raise NotImplementedError(
"AMBER currently only supports the 'full' perturbation "
"type. Please use engine='SOMD' when running multistep "
"perturbation types."
)
# If this is vacuum simulation with pmemd.cuda then
# we need to add a simulation box.
if self._is_vacuum and self._is_pmemd_cuda:
# 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))
reference_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)
reference_system = self._checkPerturbable(reference_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, 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:
# Exclude "locally enhanced sampling" executables.
if "LES" not in exe:
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)