#!/usr/bin/env python
#=============================================================================================
# MODULE DOCSTRING
#=============================================================================================
"""
Parameter handlers for the SMIRNOFF force field engine
This file contains standard parameter handlers for the SMIRNOFF force field engine.
These classes implement the object model for self-contained parameter assignment.
New pluggable handlers can be created by creating subclasses of :class:`ParameterHandler`.
.. codeauthor:: John D. Chodera <john.chodera@choderalab.org>
.. codeauthor:: David L. Mobley <dmobley@mobleylab.org>
.. codeauthor:: Peter K. Eastman <peastman@stanford.edu>
"""
#=============================================================================================
# GLOBAL IMPORTS
#=============================================================================================
from enum import Enum
import logging
from collections import OrderedDict
from simtk import openmm, unit
from openforcefield.utils import detach_units, attach_units, unit_to_string, \
extract_serialized_units_from_dict, ToolkitUnavailableException, MessageException
from openforcefield.topology import Topology, ValenceDict, ImproperDict
from openforcefield.typing.chemistry import ChemicalEnvironment
#=============================================================================================
# CONFIGURE LOGGER
#=============================================================================================
logger = logging.getLogger(__name__)
#======================================================================
# CUSTOM EXCEPTIONS
#======================================================================
class SMIRNOFFSpecError(MessageException):
"""
Exception for when data is noncompliant with the SMIRNOFF data specification.
"""
pass
class IncompatibleUnitError(MessageException):
"""
Exception for when a parameter is in the wrong units for a ParameterHandler's unit system
"""
pass
class IncompatibleParameterError(MessageException):
"""
Exception for when a set of parameters is scientifically/technically incompatible with another
"""
pass
class UnassignedValenceParameterException(Exception):
"""Exception raised when there are valence terms for which a ParameterHandler can't find parameters."""
pass
class UnassignedBondParameterException(UnassignedValenceParameterException):
"""Exception raised when there are bond terms for which a ParameterHandler can't find parameters."""
pass
class UnassignedAngleParameterException(UnassignedValenceParameterException):
"""Exception raised when there are angle terms for which a ParameterHandler can't find parameters."""
pass
class UnassignedProperTorsionParameterException(UnassignedValenceParameterException):
"""Exception raised when there are proper torsion terms for which a ParameterHandler can't find parameters."""
pass
#======================================================================
# UTILITY FUNCTIONS
#======================================================================
#======================================================================
# PARAMETER TYPE/LIST
#======================================================================
class NonbondedMethod(Enum):
"""
An enumeration of the nonbonded methods
"""
NoCutoff = 0
CutoffPeriodic = 1
CutoffNonPeriodic = 2
Ewald = 3
PME = 4
# We can't actually make this derive from dict, because it's possible for the user to change SMIRKS
# of parameters already in the list, which would cause the ParameterType object's SMIRKS and
# the dictionary key's SMIRKS to be out of sync.
[docs]class ParameterList(list):
"""
Parameter list that also supports accessing items by SMARTS string.
.. warning :: This API is experimental and subject to change.
"""
# TODO: Make this faster by caching SMARTS -> index lookup?
# TODO: Override __del__ to make sure we don't remove root atom type
# TODO: Allow retrieval by `id` as well
[docs] def __init__(self, input_parameter_list=None):
"""
Initialize a new ParameterList, optionally providing a list of ParameterType objects
to initially populate it.
Parameters
----------
input_parameter_list: list[ParameterType], default=None
A pre-existing list of ParameterType-based objects. If None, this ParameterList
will be initialized empty.
"""
super().__init__()
input_parameter_list = input_parameter_list or []
# TODO: Should a ParameterList only contain a single kind of ParameterType?
for input_parameter in input_parameter_list:
self.append(input_parameter)
[docs] def append(self, parameter):
"""
Add a ParameterType object to the end of the ParameterList
Parameters
----------
parameter : a ParameterType-derived object
"""
# TODO: Ensure that newly added parameter is the same type as existing?
super().append(parameter)
[docs] def extend(self, other):
"""
Add a ParameterList object to the end of the ParameterList
Parameters
----------
other : a ParameterList
"""
if not isinstance(other, ParameterList):
msg = 'ParameterList.extend(other) expected instance of ParameterList, ' \
'but received {} (type {}) instead'.format(other, type(other))
raise TypeError(msg)
# TODO: Check if other ParameterList contains the same ParameterTypes?
super().extend(other)
[docs] def index(self, item):
"""
Get the numerical index of a ParameterType object or SMIRKS in this ParameterList. Raises ValueError
if the item is not found.
Parameters
----------
item : ParameterType-derived object or str
The parameter or SMIRKS to look up in this ParameterList
Returns
-------
index : int
The index of the found item
"""
if isinstance(item, ParameterType):
return super().index(item)
else:
for parameter in self:
if parameter.smirks == item:
return self.index(parameter)
raise IndexError('SMIRKS {item} not found in ParameterList'.format(item=item))
[docs] def insert(self, index, parameter):
"""
Add a ParameterType object as if this were a list
Parameters
----------
index : int
The numerical position to insert the parameter at
parameter : a ParameterType-derived object
The parameter to insert
"""
# TODO: Ensure that newly added parameter is the same type as existing?
super().insert(index, parameter)
def __delitem__(self, item):
"""
Delete item by index or SMIRKS.
Parameters
----------
item : str or int
SMIRKS or numerical index of item in this ParameterList
"""
if type(item) is int:
index = item
else:
# Try to find by SMIRKS
index = self.index(item)
super().__delitem__(index)
def __getitem__(self, item):
"""
Retrieve item by index or SMIRKS
Parameters
----------
item : str or int
SMIRKS or numerical index of item in this ParameterList
"""
if type(item) is int:
index = item
elif type(item) is slice:
index = item
else:
index = self.index(item)
return super().__getitem__(index)
# TODO: Override __setitem__ and __del__ to ensure we can slice by SMIRKS as well
def __contains__(self, item):
"""Check to see if either Parameter or SMIRKS is contained in parameter list.
Parameters
----------
item : str
SMIRKS of item in this ParameterList
"""
if isinstance(item, str):
# Special case for SMIRKS strings
if item in [result.smirks for result in self]:
return True
# Fall back to traditional access
return list.__contains__(self, item)
[docs] def to_list(self, return_cosmetic_attributes=False):
"""
Render this ParameterList to a normal list, serializing each ParameterType-derived object in it to dict.
Parameters
----------
return_cosmetic_attributes : bool, optional. default = False
Whether to return non-spec attributes of each ParameterType-derived object.
Returns
-------
parameter_list : List[dict]
A serialized representation of a ParameterList, with each ParameterType it contains converted to dict.
"""
parameter_list = list()
for parameter in self:
parameter_dict = parameter.to_dict(return_cosmetic_attributes=return_cosmetic_attributes)
parameter_list.append(parameter_dict)
return parameter_list
# TODO: Rename to better reflect role as parameter base class?
[docs]class ParameterType(object):
"""
Base class for SMIRNOFF parameter types.
.. warning :: This API is experimental and subject to change.
"""
# These lists and dicts will be used for validating input and detecting cosmetic attributes.
_VALENCE_TYPE = None # ChemicalEnvironment valence type string expected by SMARTS string for this Handler
_ELEMENT_NAME = None # The string mapping to this ParameterType in a SMIRNOFF data source
_SMIRNOFF_ATTRIBS = ['smirks'] # Attributes expected per the SMIRNOFF spec.
_REQUIRE_UNITS = {} # A dict of attribs which will be checked for unit compatibility
_OPTIONAL_ATTRIBS = ['id', 'parent_id'] # Attributes in the SMIRNOFF spec that may
# be present but have no impact on performance
_INDEXED_ATTRIBS = [] # list of attribs that will have consecutive numerical suffixes starting at 1
_ATTRIBS_TO_TYPE = {} # dict of attributes that need to be cast to a type (like int or float) to be interpreted
# TODO: Can we provide some shared tools for returning settable/gettable attributes, and checking unit-bearing attributes?
[docs] def __init__(self, smirks=None, permit_cosmetic_attributes=False, **kwargs):
"""
Create a ParameterType
Parameters
----------
smirks : str
The SMIRKS match for the provided parameter type.
permit_cosmetic_attributes : bool optional. Default = False
Whether to store non-spec kwargs as "cosmetic attributes", which can be accessed and written out.
"""
import openforcefield.utils.toolkits
self._COSMETIC_ATTRIBS = [] # A list that may be populated to record the cosmetic
# attributes read from a SMIRNOFF data source
if smirks is None:
raise ValueError("'smirks' must be specified")
# TODO: Make better switch using toolkit registry
toolkit = None
if openforcefield.utils.toolkits.OPENEYE_AVAILABLE:
toolkit = 'openeye'
elif openforcefield.utils.toolkits.RDKIT_AVAILABLE:
toolkit = 'rdkit'
if toolkit is None:
raise ToolkitUnavailableException("Validating SMIRKS required either the OpenEye Toolkit or the RDKit."
" Unable to find either.")
ChemicalEnvironment.validate(
smirks, ensure_valence_type=self._VALENCE_TYPE, toolkit=toolkit)
self._smirks = smirks
def _assert_quantity_is_compatible(quantity_name, quantity, unit_to_check):
"""
Checks whether a simtk.unit.Quantity is compatible with another unit.
Parameters
----------
quantity_name : string
quantity : A simtk.unit.Quantity
unit_to_check : A simtk.unit.Unit
"""
if not quantity.unit.is_compatible(unit_to_check):
msg = "{} constructor received kwarg {} with value {}, " \
"which is incompatible with expected unit {}".format(self.__class__,
quantity_name,
val,
unit_to_check)
raise SMIRNOFFSpecError(msg)
# First look for indexed attribs, removing them from kwargs as they're found
for unidx_key in self._INDEXED_ATTRIBS:
# Start by generating a string with the key + an index
index = 1
idx_key = unidx_key+str(index)
# If the indexed key is present in the kwargs, set the attrib data type to be a list
if idx_key in kwargs:
setattr(self, unidx_key, list())
# Iterate through increasing values on the index, appending them as they are found
while idx_key in kwargs:
val = kwargs[idx_key]
# If the indexed keys require units, ensure they are compatible
if unidx_key in self._REQUIRE_UNITS:
_assert_quantity_is_compatible(idx_key, val, self._REQUIRE_UNITS[unidx_key])
# If the indexed keys need to be cast to a type, do that here
if unidx_key in self._ATTRIBS_TO_TYPE:
type_to_cast = self._ATTRIBS_TO_TYPE[unidx_key]
val = type_to_cast(val)
# Finally, append this value to the attribute and remove the key from kwargs
getattr(self, unidx_key).append(val)
del kwargs[idx_key]
index += 1
idx_key = unidx_key + str(index)
# Iterate through the remaining, non-indexed kwargs,
# doing validation and setting this ParameterType's attributes
for key, val in kwargs.items():
if key in self._REQUIRE_UNITS:
# TODO: Add dynamic property-getter/setter for each thing in self._REQUIRE_UNITS
_assert_quantity_is_compatible(key, val, self._REQUIRE_UNITS[key])
if key in self._ATTRIBS_TO_TYPE:
type_to_cast = self._ATTRIBS_TO_TYPE[key]
val = type_to_cast(val)
# Iterate through
# TODO: Decide on a scheme for prepending underscores to attributes
#attr_name = '_' + key
if key in self._SMIRNOFF_ATTRIBS:
setattr(self, key, val)
# If it's an optional attrib,
elif key in self._OPTIONAL_ATTRIBS:
setattr(self, key, val)
# Handle all unknown kwargs as cosmetic so we can write them back out
elif permit_cosmetic_attributes:
self._COSMETIC_ATTRIBS.append(key)
setattr(self, key, val)
else:
raise SMIRNOFFSpecError("Incompatible kwarg {} passed to {} constructor. If this is "
"a desired cosmetic attribute, consider setting "
"'permit_cosmetic_attributes=True'".format({key: val}, self.__class__))
@property
def smirks(self):
return self._smirks
@smirks.setter
def smirks(self, smirks):
# Validate the SMIRKS string to ensure it matches the expected parameter type,
# raising an exception if it is invalid or doesn't tag a valid set of atoms
# TODO: Add check to make sure we can't make tree non-hierarchical
# This would require parameter type knows which ParameterList it belongs to
ChemicalEnvironment.validate(
smirks, ensure_valence_type=self._VALENCE_TYPE)
self._smirks = smirks
[docs] def to_dict(self, return_cosmetic_attributes=False):
"""
Convert this ParameterType-derived object to dict. A unit-bearing attribute ('X') will be converted to two dict
entries, one (['X'] containing the unitless value, and another (['X_unit']) containing a string representation
of its unit.
Parameters
----------
return_cosmetic_attributes : bool, optional. default = False
Whether to return non-spec attributes of this ParameterType
Returns
-------
smirnoff_dict : dict
The SMIRNOFF-compliant dict representation of this ParameterType-derived object.
output_units : dict[str: simtk.unit.Unit]
A mapping from each simtk.unit.Quanitity-valued ParameterType attribute
to the unit it was converted to during serialization.
"""
# Make a list of all attribs that should be included in the
# returned dict (call list() to make a copy)
attribs_to_return = list(self._SMIRNOFF_ATTRIBS)
attribs_to_return += [opt_attrib for opt_attrib in self._OPTIONAL_ATTRIBS if hasattr(self, opt_attrib)]
if return_cosmetic_attributes:
attribs_to_return += self._COSMETIC_ATTRIBS
# Start populating a dict of the attribs
smirnoff_dict = OrderedDict()
# If attribs_to_return is ordered here, that will effectively be an informal output ordering
for attrib_name in attribs_to_return:
attrib_value = self.__getattribute__(attrib_name)
if type(attrib_value) is list:
for idx, val in enumerate(attrib_value):
smirnoff_dict[attrib_name + str(idx+1)] = val
else:
smirnoff_dict[attrib_name] = attrib_value
return smirnoff_dict
def __repr__(self):
ret_str = '<{} with '.format(self.__class__.__name__)
for attr, val in self.to_dict().items():
ret_str += f'{attr}: {val} '
ret_str += '>'
return ret_str
#======================================================================
# PARAMETER HANDLERS
#
# The following classes are Handlers that know how to create Force
# subclasses and add them to a System that is being created. Each Handler
# class must define three methods:
# 1) a constructor which takes as input hierarchical dictionaries of data
# conformant to the SMIRNOFF spec;
# 2) a create_force() method that constructs the Force object and adds it
# to the System; and
# 3) a labelForce() method that provides access to which terms are applied
# to which atoms in specified mols.
#======================================================================
# TODO: Should we have a parameter handler registry?
[docs]class ParameterHandler(object):
"""Base class for parameter handlers.
Parameter handlers are configured with some global parameters for a given section. They may also contain a
:class:`ParameterList` populated with :class:`ParameterType` objects if they are responsile for assigning
SMIRKS-based parameters.
.. warning
Parameter handler objects can only belong to a single :class:`ForceField` object.
If you need to create a copy to attach to a different :class:`ForceField` object, use ``create_copy()``.
.. warning :: This API is experimental and subject to change.
"""
_TAGNAME = None # str of section type handled by this ParameterHandler (XML element name for SMIRNOFF XML representation)
_INFOTYPE = None # container class with type information that will be stored in self._parameters
_OPENMMTYPE = None # OpenMM Force class (or None if no equivalent)
_DEPENDENCIES = None # list of ParameterHandler classes that must precede this, or None
_REQUIRED_SPEC_ATTRIBS = [] # list of kwargs that must be present during handler initialization
_DEFAULT_SPEC_ATTRIBS = {} # dict of tag-level attributes and their default values
_OPTIONAL_SPEC_ATTRIBS = [] # list of non-required attributes that can be defined on initialization
_INDEXED_ATTRIBS = [] # list of parameter attribs that will have consecutive numerical suffixes starting at 1
_REQUIRE_UNITS = {} # dict of {header attrib : unit } for input checking
_ATTRIBS_TO_TYPE = {} # dict of attribs that must be cast to a specific type to be interpreted correctly
_KWARGS = [] # Kwargs to catch when create_force is called
_SMIRNOFF_VERSION_INTRODUCED = 0.0 # the earliest version of SMIRNOFF spec that supports this ParameterHandler
_SMIRNOFF_VERSION_DEPRECATED = None # if deprecated, the first SMIRNOFF version number it is no longer used
[docs] def __init__(self, permit_cosmetic_attributes=False, **kwargs):
"""
Initialize a ParameterHandler, optionally with a list of parameters and other kwargs.
Parameters
----------
permit_cosmetic_attributes : bool
Whether to accept non-spec kwargs
**kwargs : dict
The dict representation of the SMIRNOFF data source
"""
self._COSMETIC_ATTRIBS = [] # list of cosmetic header attributes to remember and optionally write out
# Ensure that all required attribs are present
for reqd_attrib in self._REQUIRED_SPEC_ATTRIBS:
if not reqd_attrib in kwargs:
msg = "{} requires {} as a parameter during initialization, however this is not " \
"provided. Defined kwargs are {}".format(self.__class__,
reqd_attrib,
list(kwargs.keys()))
raise SMIRNOFFSpecError(msg)
# list of ParameterType objects (also behaves like an OrderedDict where keys are SMARTS)
self._parameters = ParameterList()
# Handle all the unknown kwargs as cosmetic so we can write them back out
allowed_header_attribs = self._REQUIRED_SPEC_ATTRIBS + \
list(self._DEFAULT_SPEC_ATTRIBS.keys()) + \
self._OPTIONAL_SPEC_ATTRIBS
# Check for attribs that need to be casted to specific types
for attrib, type_to_cast in self._ATTRIBS_TO_TYPE.items():
if attrib in kwargs:
kwargs[attrib] = type_to_cast(kwargs[attrib])
# Check for indexed attribs
for attrib_basename in self._INDEXED_ATTRIBS:
attrib_unit_key = attrib_basename + '_unit'
index = 1
attrib_w_index = '{}{}'.format(attrib_basename, index)
while attrib_w_index in kwargs:
# As long as we keep finding higher-indexed entries for
# this attrib, add them to the expected arguments
allowed_header_attribs.append(attrib_w_index)
# If there's a unit for this attrib, copy unit entries for each index instance
if attrib_unit_key in kwargs:
kwargs[attrib_w_index+'_unit'] = kwargs[attrib_unit_key]
# Attach units to the handler kwargs, if applicable
unitless_kwargs, attached_units = extract_serialized_units_from_dict(kwargs)
smirnoff_data = attach_units(unitless_kwargs, attached_units)
# Add default values to smirnoff_data if they're not already there
for default_key, default_val in self._DEFAULT_SPEC_ATTRIBS.items():
if not (default_key in kwargs):
smirnoff_data[default_key] = default_val
# Perform unit compatibility checks
for key, val in smirnoff_data.items():
if key in self._REQUIRE_UNITS:
# TODO: Logic for indexed attributes
if not val.unit.is_compatible(self._REQUIRE_UNITS[key]):
msg = "{} constructor received kwarg {} with value {}, " \
"which is incompatible with expected unit {}".format(self.__class__,
key,
val,
self._REQUIRE_UNITS[key])
raise SMIRNOFFSpecError(msg)
element_name = None
if self._INFOTYPE is not None:
element_name = self._INFOTYPE._ELEMENT_NAME
for key, val in smirnoff_data.items():
# If we're reading the parameter list, iterate through and attach units to
# each parameter_dict, then use it to initialize a ParameterType
if key == element_name:
# If there are multiple parameters, this will be a list. If there's just one, make it a list
if not(isinstance(val, list)):
val = [val]
for unitless_param_dict in val:
param_dict = attach_units(unitless_param_dict, attached_units)
new_parameter = self._INFOTYPE(**param_dict,
permit_cosmetic_attributes=permit_cosmetic_attributes)
self._parameters.append(new_parameter)
elif key in allowed_header_attribs:
attr_name = '_' + key
# TODO: create @property.setter here if attrib requires unit
setattr(self, attr_name, val)
elif permit_cosmetic_attributes:
self._COSMETIC_ATTRIBS.append(key)
attr_name = '_' + key
setattr(self, attr_name, val)
else:
raise SMIRNOFFSpecError("Incompatible kwarg {} passed to {} constructor. If this is "
"a desired cosmetic attribute, consider setting "
"'permit_cosmetic_attributes=True'".format(key, self.__class__))
@property
def parameters(self):
"""The ParameterList that holds this ParameterHandler's parameter objects"""
return self._parameters
# TODO: Do we need to return these, or can we handle this internally
@property
def known_kwargs(self):
"""List of kwargs that can be parsed by the function.
"""
# TODO: Should we use introspection to inspect the function signature instead?
return set(self._KWARGS)
#@classmethod
[docs] def check_parameter_compatibility(self, parameter_kwargs):
"""
Check to make sure that the fields requiring defined units are compatible with the required units for the
Parameters handled by this ParameterHandler
Parameters
----------
parameter_kwargs: dict
The dict that will be used to construct the ParameterType
Raises
------
Raises a ValueError if the parameters are incompatible.
"""
for key in parameter_kwargs:
if key in self._REQUIRE_UNITS:
reqd_unit = self._REQUIRE_UNITS[key]
#if arg in cls._REQUIRE_UNITS:
# raise Exception(cls)
# reqd_unit = cls._REQUIRE_UNITS[arg]
val = parameter_kwargs[key]
if not (reqd_unit.is_compatible(val.unit)):
raise IncompatibleUnitError(
"Input unit {} is not compatible with ParameterHandler unit {}"
.format(val.unit, reqd_unit))
[docs] def check_handler_compatibility(self, handler_kwargs):
"""
Checks if a set of kwargs used to create a ParameterHandler are compatible with this ParameterHandler. This is
called if a second handler is attempted to be initialized for the same tag.
Parameters
----------
handler_kwargs : dict
The kwargs that would be used to construct
Raises
------
IncompatibleParameterError if handler_kwargs are incompatible with existing parameters.
"""
pass
# TODO: Can we ensure SMIRKS and other parameters remain valid after manipulation?
[docs] def add_parameter(self, parameter_kwargs):
"""Add a parameter to the forcefield, ensuring all parameters are valid.
Parameters
----------
parameter_kwargs : dict
The kwargs to pass to the ParameterHandler.INFOTYPE (a ParameterType) constructor
"""
# TODO: Do we need to check for incompatibility with existing parameters?
# Perform unit compatibility checks
self.check_parameter_compatibility(parameter_kwargs)
# Check for correct SMIRKS valence
new_parameter = self._INFOTYPE(**parameter_kwargs)
self._parameters.append(new_parameter)
[docs] def get_parameter(self, parameter_attrs):
"""
Return the parameters in this ParameterHandler that match the parameter_attrs argument
Parameters
----------
parameter_attrs : dict of {attr: value}
The attrs mapped to desired values (for example {"smirks": "[*:1]~[#16:2]=,:[#6:3]~[*:4]", "id": "t105"} )
Returns
-------
list of ParameterType-derived objects
A list of matching ParameterType-derived objects
"""
# TODO: This is a necessary API point for Lee-Ping's ForceBalance
pass
[docs] def find_matches(self, entity):
"""Find the elements of the topology/molecule matched by a parameter type.
Parameters
----------
entity : openforcefield.topology.Topology or openforcefield.topology.Molecule
Topology or molecule to search.
Returns
---------
matches : ValenceDict[Tuple[int], ParameterType]
``matches[particle_indices]`` is the ``ParameterType`` object
matching the tuple of particle indices in ``entity``.
"""
return self._find_matches(entity)
def _find_matches(self, entity, transformed_dict_cls=ValenceDict):
"""Implement find_matches() and allow using a difference valence dictionary."""
from openforcefield.topology import FrozenMolecule
logger.debug('Finding matches for {}'.format(self.__class__.__name__))
matches = transformed_dict_cls()
for parameter_type in self._parameters:
matches_for_this_type = {}
for atoms in entity.chemical_environment_matches(parameter_type.smirks):
# Collect the atom indices matching the entity.
if isinstance(entity, Topology):
atom_indices = tuple([atom.topology_particle_index for atom in atoms])
elif isinstance(entity, FrozenMolecule):
atom_indices = tuple([atom.molecule_particle_index for atom in atoms])
else:
raise ValueError('Unknown entity type {}'.format(entity.__class__))
# Update the matches for this parameter type.
matches_for_this_type[atom_indices] = parameter_type
# Update matches of all parameter types.
matches.update(matches_for_this_type)
logger.debug('{:64} : {:8} matches'.format(
parameter_type.smirks, len(matches_for_this_type)))
logger.debug('{} matches identified'.format(len(matches)))
return matches
[docs] def assign_parameters(self, topology, system):
"""Assign parameters for the given Topology to the specified System object.
Parameters
----------
topology : openforcefield.topology.Topology
The Topology for which parameters are to be assigned.
Either a new Force will be created or parameters will be appended to an existing Force.
system : simtk.openmm.System
The OpenMM System object to add the Force (or append new parameters) to.
"""
pass
[docs] def postprocess_system(self, topology, system, **kwargs):
"""Allow the force to perform a a final post-processing pass on the System following parameter assignment, if needed.
Parameters
----------
topology : openforcefield.topology.Topology
The Topology for which parameters are to be assigned.
Either a new Force will be created or parameters will be appended to an existing Force.
system : simtk.openmm.System
The OpenMM System object to add the Force (or append new parameters) to.
"""
pass
[docs] def to_dict(self, output_units=None, return_cosmetic_attributes=False):
"""
Convert this ParameterHandler to an OrderedDict, compliant with the SMIRNOFF data spec.
Parameters
----------
output_units : dict[str : simtk.unit.Unit], optional. Default = None
A mapping from the ParameterType attribute name to the output unit its value should be converted to.
return_cosmetic_attributes : bool, optional. Default = False.
Whether to return non-spec parameter and header attributes in this ParameterHandler.
Returns
-------
smirnoff_data : OrderedDict
SMIRNOFF-spec compliant representation of this ParameterHandler and its internal ParameterList.
"""
smirnoff_data = OrderedDict()
# Set default output units to those from the last parameter added to the ParameterList
if (output_units is None):
output_units = dict()
# TODO: What if self._parameters is an empty list?
if len(self._parameters) > 0:
_, last_added_output_units = detach_units(self._parameters[-1].to_dict())
# Overwrite key_value pairs in last_added_output_units with those specified by user in output_units
last_added_output_units.update(output_units)
output_units = last_added_output_units
# Populate parameter list
parameter_list = self._parameters.to_list(return_cosmetic_attributes=return_cosmetic_attributes)
unitless_parameter_list = list()
# Detach units into a separate dict.
for parameter_dict in parameter_list:
unitless_parameter_dict, attached_units = detach_units(parameter_dict, output_units=output_units)
unitless_parameter_list.append(unitless_parameter_dict)
output_units.update(attached_units)
# Collapse down indexed attribute units
# (eg. {'k1_unit': angstrom, 'k2_unit': angstrom} --> {'k_unit': angstrom})
for attrib_key in self._INDEXED_ATTRIBS:
index = 1
# Store a variable that is 'k1_unit'
idxed_attrib_unit_key = attrib_key + str(index) + '_unit'
# See if 'k1_unit' is in output_units
if idxed_attrib_unit_key in output_units:
# If so, define 'k_unit' and add it to the output_units dict
attrib_unit_key = attrib_key + '_unit'
output_units[attrib_unit_key] = output_units[idxed_attrib_unit_key]
# Increment the 'kN_unit' value, checking that each is the same as the
# 'k1_unit' value, and deleting them from output_units
while idxed_attrib_unit_key in output_units:
# Ensure that no different units are defined for higher indexes of this attrib
assert output_units[attrib_unit_key] == output_units[idxed_attrib_unit_key]
del output_units[idxed_attrib_unit_key]
index += 1
idxed_attrib_unit_key = attrib_key + str(index) + '_unit'
# NOTE: This assumes that a ParameterHandler will have just one homogenous ParameterList under it
if self._INFOTYPE is not None:
smirnoff_data[self._INFOTYPE._ELEMENT_NAME] = unitless_parameter_list
# Collect the names of handler attributes to return
header_attribs_to_return = self._REQUIRED_SPEC_ATTRIBS + list(self._DEFAULT_SPEC_ATTRIBS.keys())
# Check whether the optional attribs are defined, and add them if so
for key in self._OPTIONAL_SPEC_ATTRIBS:
attr_key = '_' + key
if hasattr(self, attr_key):
header_attribs_to_return.append(key)
# Add the cosmetic attributes if requested
if return_cosmetic_attributes:
header_attribs_to_return += self._COSMETIC_ATTRIBS
# Go through the attribs of this ParameterHandler and collect the appropriate values to return
header_attribute_dict = {}
for header_attribute in header_attribs_to_return:
value = getattr(self, '_' + header_attribute)
header_attribute_dict[header_attribute] = value
# Detach all units from the header attribs
unitless_header_attribute_dict, attached_header_units = detach_units(header_attribute_dict)
# Convert all header attrib units (eg. {'length_unit': simtk.unit.angstrom}) to strings (eg.
# {'length_unit': 'angstrom'}) and add them to the header attribute dict
# TODO: Should we check for collisions between parameter "_unit" keys and header "_unit" keys?
output_units.update(attached_header_units)
for key, value in output_units.items():
value_str = unit_to_string(value)
# Made a note to add this to the smirnoff spec
output_units[key] = value_str
# Reattach the attached units here
smirnoff_data.update(unitless_header_attribute_dict)
smirnoff_data.update(output_units)
return smirnoff_data
# -------------------------------
# Utilities for children classes.
# -------------------------------
@classmethod
def _check_all_valence_terms_assigned(cls, assigned_terms, valence_terms,
exception_cls=UnassignedValenceParameterException):
"""Check that all valence terms have been assigned and print a user-friendly error message.
Parameters
----------
assigned_terms : ValenceDict
Atom index tuples defining added valence terms.
valence_terms : Iterable[TopologyAtom] or Iterable[Iterable[TopologyAtom]]
Atom or atom tuples defining topological valence terms.
exception_cls : UnassignedValenceParameterException
A specific exception class to raise to allow catching only specific
types of errors.
"""
# Convert the valence term to a valence dictionary to make sure
# the order of atom indices doesn't matter for comparison.
valence_terms_dict = assigned_terms.__class__()
for atoms in valence_terms:
try:
# valence_terms is a list of TopologyAtom tuples.
atom_indices = (a.topology_particle_index for a in atoms)
except TypeError:
# valence_terms is a list of TopologyAtom.
atom_indices = (atoms.topology_particle_index,)
valence_terms_dict[atom_indices] = atoms
# Check that both valence dictionaries have the same keys (i.e. terms).
assigned_terms_set = set(assigned_terms.keys())
valence_terms_set = set(valence_terms_dict.keys())
unassigned_terms = valence_terms_set.difference(assigned_terms_set)
not_found_terms = assigned_terms_set.difference(valence_terms_set)
# Raise an error if there are unassigned terms.
err_msg = ""
if len(unassigned_terms) > 0:
unassigned_str = '\n- '.join([str(x) for x in unassigned_terms])
err_msg += ("{parameter_handler} was not able to find parameters for the following valence terms:\n"
"- {unassigned_str}").format(parameter_handler=cls.__name__,
unassigned_str=unassigned_str)
if len(not_found_terms) > 0:
if err_msg != "":
err_msg += '\n'
not_found_str = '\n- '.join([str(x) for x in not_found_terms])
err_msg += ("{parameter_handler} assigned terms that were not found in the topology:\n"
"- {not_found_str}").format(parameter_handler=cls.__name__,
not_found_str=not_found_str)
if err_msg != "":
err_msg += '\n'
raise exception_cls(err_msg)
#=============================================================================================
class ConstraintHandler(ParameterHandler):
"""Handle SMIRNOFF ``<Constraints>`` tags
``ConstraintHandler`` must be applied before ``BondHandler`` and ``AngleHandler``,
since those classes add constraints for which equilibrium geometries are needed from those tags.
.. warning :: This API is experimental and subject to change.
"""
class ConstraintType(ParameterType):
"""A SMIRNOFF constraint type
.. warning :: This API is experimental and subject to change.
"""
_VALENCE_TYPE = 'Bond'
_ELEMENT_NAME = 'Constraint'
_SMIRNOFF_ATTRIBS = ['smirks'] # Attributes expected per the SMIRNOFF spec.
_OPTIONAL_ATTRIBS = ['distance', 'id', 'parent_id']
_REQUIRE_UNITS = {'distance': unit.angstrom}
def __init__(self, **kwargs):
super().__init__(**kwargs)
# TODO: Re-implement ability to set 'distance = True'
# if 'distance' in node.attrib:
# self.distance = _extract_quantity_from_xml_element(
# node, parent, 'distance'
# ) # Constraint with specified distance will be added by ConstraintHandler
# else:
# self.distance = True # Constraint to equilibrium bond length will be added by HarmonicBondHandler
_TAGNAME = 'Constraints'
_INFOTYPE = ConstraintType
_OPENMMTYPE = None # don't create a corresponding OpenMM Force class
def __init__(self, **kwargs):
super().__init__(**kwargs)
def create_force(self, system, topology, **kwargs):
constraints = self.find_matches(topology)
for (atoms, constraint) in constraints.items():
# Update constrained atom pairs in topology
#topology.add_constraint(*atoms, constraint.distance)
# If a distance is specified (constraint.distance != True), add the constraint here.
# Otherwise, the equilibrium bond length will be used to constrain the atoms in HarmonicBondHandler
if hasattr(constraint, 'distance'):# is not True:
system.addConstraint(*atoms, constraint.distance)
topology.add_constraint(*atoms, constraint.distance)
else:
topology.add_constraint(*atoms, True)
#=============================================================================================
[docs]class BondHandler(ParameterHandler):
"""Handle SMIRNOFF ``<Bonds>`` tags
.. warning :: This API is experimental and subject to change.
"""
[docs] class BondType(ParameterType):
"""A SMIRNOFF bond type
.. warning :: This API is experimental and subject to change.
"""
_VALENCE_TYPE = 'Bond' # ChemicalEnvironment valence type string expected by SMARTS string for this Handler
_ELEMENT_NAME = 'Bond'
_SMIRNOFF_ATTRIBS = ['smirks', 'length', 'k'] # Attributes expected per the SMIRNOFF spec.
_REQUIRE_UNITS = {'length' : unit.angstrom,
'k' : unit.kilocalorie_per_mole / unit.angstrom**2}
_INDEXED_ATTRIBS = ['length', 'k'] # May be indexed (by integer bond order) if fractional bond orders are used
def __init__(self, **kwargs):
super().__init__(**kwargs) # Base class handles ``smirks`` and ``id`` fields
_TAGNAME = 'Bonds' # SMIRNOFF tag name to process
_INFOTYPE = BondType # class to hold force type info
_OPENMMTYPE = openmm.HarmonicBondForce # OpenMM force class to create
_DEPENDENCIES = [ConstraintHandler] # ConstraintHandler must be executed first
_DEFAULT_SPEC_ATTRIBS = {'potential': 'harmonic',
'fractional_bondorder_method': None,
'fractional_bondorder_interpolation': 'linear'}
_INDEXED_ATTRIBS = ['k'] # May be indexed (by integer bond order) if fractional bond orders are used
[docs] def __init__(self, **kwargs):
# TODO: Do we want a docstring here? If not, check that docstring get inherited from ParameterHandler.
super().__init__(**kwargs)
def create_force(self, system, topology, **kwargs):
# Create or retrieve existing OpenMM Force object
# TODO: The commented line below should replace the system.getForce search
#force = super(BondHandler, self).create_force(system, topology, **kwargs)
existing = [system.getForce(i) for i in range(system.getNumForces())]
existing = [f for f in existing if type(f) == self._OPENMMTYPE]
if len(existing) == 0:
force = self._OPENMMTYPE()
system.addForce(force)
else:
force = existing[0]
# Add all bonds to the system.
bonds = self.find_matches(topology)
skipped_constrained_bonds = 0 # keep track of how many bonds were constrained (and hence skipped)
for (atoms, bond_params) in bonds.items():
# Get corresponding particle indices in Topology
#particle_indices = tuple([ atom.particle_index for atom in atoms ])
# Ensure atoms are actually bonded correct pattern in Topology
topology.assert_bonded(atoms[0], atoms[1])
# Compute equilibrium bond length and spring constant.
topology_bond = topology.get_bond_between(atoms[0], atoms[1])
if topology_bond.bond.fractional_bond_order is None:
[k, length] = [bond_params.k, bond_params.length]
else:
# Interpolate using fractional bond orders
# TODO: Do we really want to allow per-bond specification of interpolation schemes?
order = topology_bond.bond.fractional_bond_order
if self.fractional_bondorder_interpolation == 'interpolate-linear':
k = bond_params.k[0] + (bond_params.k[1] - bond_params.k[0]) * (order - 1.)
length = bond_params.length[0] + (
bond_params.length[1] - bond_params.length[0]) * (order - 1.)
else:
raise Exception(
"Partial bondorder treatment {} is not implemented.".
format(self.fractional_bondorder_method))
# Handle constraints.
if topology.is_constrained(*atoms):
# Atom pair is constrained; we don't need to add a bond term.
skipped_constrained_bonds += 1
# Check if we need to add the constraint here to the equilibrium bond length.
if topology.is_constrained(*atoms) is True:
# Mark that we have now assigned a specific constraint distance to this constraint.
topology.add_constraint(*atoms, length)
# Add the constraint to the System.
system.addConstraint(*atoms, length)
#system.addConstraint(*particle_indices, length)
continue
# Add harmonic bond to HarmonicBondForce
force.addBond(*atoms, length, k)
logger.info('{} bonds added ({} skipped due to constraints)'.format(
len(bonds) - skipped_constrained_bonds, skipped_constrained_bonds))
# Check that no topological bonds are missing force parameters.
valence_terms = [list(b.atoms) for b in topology.topology_bonds]
self._check_all_valence_terms_assigned(assigned_terms=bonds, valence_terms=valence_terms,
exception_cls=UnassignedBondParameterException)
#=============================================================================================
[docs]class AngleHandler(ParameterHandler):
"""Handle SMIRNOFF ``<AngleForce>`` tags
.. warning :: This API is experimental and subject to change.
"""
[docs] class AngleType(ParameterType):
"""A SMIRNOFF angle type.
.. warning :: This API is experimental and subject to change.
"""
_VALENCE_TYPE = 'Angle' # ChemicalEnvironment valence type string expected by SMARTS string for this Handler
_ELEMENT_NAME = 'Angle'
_SMIRNOFF_ATTRIBS = ['smirks', 'angle', 'k'] # Attributes expected per the SMIRNOFF spec.
_REQUIRE_UNITS = {'angle': unit.degree,
'k': unit.kilocalorie_per_mole / unit.degree**2}
def __init__(self, **kwargs):
super().__init__(**kwargs) # base class handles ``smirks`` and ``id`` fields
_TAGNAME = 'Angles' # SMIRNOFF tag name to process
_INFOTYPE = AngleType # class to hold force type info
_OPENMMTYPE = openmm.HarmonicAngleForce # OpenMM force class to create
_DEFAULT_SPEC_ATTRIBS = {'potential': 'harmonic'}
[docs] def __init__(self, **kwargs):
super().__init__(**kwargs)
def create_force(self, system, topology, **kwargs):
#force = super(AngleHandler, self).create_force(system, topology, **kwargs)
existing = [system.getForce(i) for i in range(system.getNumForces())]
existing = [f for f in existing if type(f) == self._OPENMMTYPE]
if len(existing) == 0:
force = self._OPENMMTYPE()
system.addForce(force)
else:
force = existing[0]
# Add all angles to the system.
angles = self.find_matches(topology)
skipped_constrained_angles = 0 # keep track of how many angles were constrained (and hence skipped)
for (atoms, angle) in angles.items():
# Ensure atoms are actually bonded correct pattern in Topology
for (i, j) in [(0, 1), (1, 2)]:
topology.assert_bonded(atoms[i], atoms[j])
if topology.is_constrained(
atoms[0], atoms[1]) and topology.is_constrained(
atoms[1], atoms[2]) and topology.is_constrained(
atoms[0], atoms[2]):
# Angle is constrained; we don't need to add an angle term.
skipped_constrained_angles += 1
continue
force.addAngle(*atoms, angle.angle, angle.k)
logger.info('{} angles added ({} skipped due to constraints)'.format(
len(angles) - skipped_constrained_angles,
skipped_constrained_angles))
# Check that no topological angles are missing force parameters
self._check_all_valence_terms_assigned(assigned_terms=angles,
valence_terms=list(topology.angles),
exception_cls=UnassignedAngleParameterException)
#=============================================================================================
[docs]class ProperTorsionHandler(ParameterHandler):
"""Handle SMIRNOFF ``<ProperTorsionForce>`` tags
.. warning :: This API is experimental and subject to change.
"""
[docs] class ProperTorsionType(ParameterType):
"""A SMIRNOFF torsion type for proper torsions.
.. warning :: This API is experimental and subject to change.
"""
_VALENCE_TYPE = 'ProperTorsion'
_ELEMENT_NAME = 'Proper'
_SMIRNOFF_ATTRIBS = ['smirks', 'periodicity', 'phase', 'k'] # Attributes expected per the SMIRNOFF spec.
_REQUIRE_UNITS = {'k': unit.kilocalorie_per_mole,
'phase': unit.degree}
_OPTIONAL_ATTRIBS = ['id', 'parent_id', 'idivf']
_INDEXED_ATTRIBS = ['k', 'phase', 'periodicity', 'idivf']
_ATTRIBS_TO_TYPE = {'periodicity': int,
'idivf': float}
def __init__(self, **kwargs):
super().__init__(**kwargs) # base class handles ``smirks`` and ``id`` fields
_TAGNAME = 'ProperTorsions' # SMIRNOFF tag name to process
_INFOTYPE = ProperTorsionType # info type to store
_OPENMMTYPE = openmm.PeriodicTorsionForce # OpenMM force class to create
_DEFAULT_SPEC_ATTRIBS = {'potential': 'charmm',
'default_idivf': 'auto'}
_INDEXED_ATTRIBS = ['k', 'phase', 'periodicity', 'idivf']
[docs] def __init__(self, **kwargs):
# NOTE: We do not want to overwrite idivf values here! If they're missing from the ParameterType
# dictionary, that means they should be set to defualt _AT SYSTEM CREATION TIME_. The user may
# change that default to a different value than it is now. The solution here will be to leave
# those idivfX values uninitialized and deal with it during system creation
super().__init__(**kwargs)
def create_force(self, system, topology, **kwargs):
#force = super(ProperTorsionHandler, self).create_force(system, topology, **kwargs)
existing = [system.getForce(i) for i in range(system.getNumForces())]
existing = [f for f in existing if type(f) == self._OPENMMTYPE]
if len(existing) == 0:
force = self._OPENMMTYPE()
system.addForce(force)
else:
force = existing[0]
# Add all proper torsions to the system.
torsions = self.find_matches(topology)
for (atom_indices, torsion) in torsions.items():
# Ensure atoms are actually bonded correct pattern in Topology
for (i, j) in [(0, 1), (1, 2), (2, 3)]:
topology.assert_bonded(atom_indices[i], atom_indices[j])
for (periodicity, phase, k, idivf) in zip(torsion.periodicity,
torsion.phase, torsion.k, torsion.idivf):
if idivf == 'auto':
# TODO: Implement correct "auto" behavior
raise NotImplementedError("The OpenForceField toolkit hasn't implemented "
"support for the torsion `idivf` value of 'auto'")
force.addTorsion(atom_indices[0], atom_indices[1],
atom_indices[2], atom_indices[3], periodicity,
phase, k/idivf)
logger.info('{} torsions added'.format(len(torsions)))
# Check that no topological torsions are missing force parameters
self._check_all_valence_terms_assigned(assigned_terms=torsions,
valence_terms=list(topology.propers),
exception_cls=UnassignedProperTorsionParameterException)
[docs]class ImproperTorsionHandler(ParameterHandler):
"""Handle SMIRNOFF ``<ImproperTorsionForce>`` tags
.. warning :: This API is experimental and subject to change.
"""
[docs] class ImproperTorsionType(ParameterType):
"""A SMIRNOFF torsion type for improper torsions.
.. warning :: This API is experimental and subject to change.
"""
_VALENCE_TYPE = 'ImproperTorsion'
_ELEMENT_NAME = 'Improper'
_SMIRNOFF_ATTRIBS = ['smirks', 'periodicity', 'phase', 'k'] # Attributes expected per the SMIRNOFF spec.
_REQUIRE_UNITS = {'k': unit.kilocalorie_per_mole,
'phase': unit.degree}
_OPTIONAL_ATTRIBS = ['id', 'parent_id', 'idivf']
_INDEXED_ATTRIBS = ['k', 'phase', 'periodicity', 'idivf']
_ATTRIBS_TO_TYPE = {'periodicity': int,
'idivf': float}
def __init__(self, **kwargs):
super().__init__( **kwargs)
_TAGNAME = 'ImproperTorsions' # SMIRNOFF tag name to process
_INFOTYPE = ImproperTorsionType # info type to store
_OPENMMTYPE = openmm.PeriodicTorsionForce # OpenMM force class to create
_OPTIONAL_SPEC_ATTRIBS = ['potential', 'default_idivf']
_HANDLER_DEFAULTS = {'potential': 'charmm',
'default_idivf': 'auto'}
_INDEXED_ATTRIBS = ['k', 'phase', 'periodicity', 'idivf']
[docs] def __init__(self, **kwargs):
super().__init__(**kwargs)
[docs] def find_matches(self, entity):
"""Find the improper torsions in the topology/molecule matched by a parameter type.
Parameters
----------
entity : openforcefield.topology.Topology or openforcefield.topology.Molecule
Topology or molecule to search.
Returns
---------
matches : ImproperDict[Tuple[int], ParameterType]
``matches[atom_indices]`` is the ``ParameterType`` object
matching the 4-tuple of atom indices in ``entity``.
"""
return self._find_matches(entity, transformed_dict_cls=ImproperDict)
def create_force(self, system, topology, **kwargs):
#force = super(ImproperTorsionHandler, self).create_force(system, topology, **kwargs)
#force = super().create_force(system, topology, **kwargs)
existing = [system.getForce(i) for i in range(system.getNumForces())]
existing = [
f for f in existing if type(f) == openmm.PeriodicTorsionForce
]
if len(existing) == 0:
force = openmm.PeriodicTorsionForce()
system.addForce(force)
else:
force = existing[0]
# Add all improper torsions to the system
impropers = self.find_matches(topology)
for (atom_indices, improper) in impropers.items():
# Ensure atoms are actually bonded correct pattern in Topology
# For impropers, central atom is atom 1
for (i, j) in [(0, 1), (1, 2), (1, 3)]:
topology.assert_bonded(atom_indices[i], atom_indices[j])
# TODO: This is a lazy hack. idivf should be set according to the ParameterHandler's default_idivf attrib
if not hasattr(improper, 'idivf'):
improper.idivf = [3 for item in improper.k]
# Impropers are applied in three paths around the trefoil having the same handedness
for (improper_periodicity, improper_phase, improper_k, improper_idivf) in zip(improper.periodicity,
improper.phase, improper.k, improper.idivf):
# TODO: Implement correct "auto" behavior
if improper_idivf == 'auto':
improper_idivf = 3
logger.warning("The OpenForceField toolkit hasn't implemented "
"support for the torsion `idivf` value of 'auto'."
"Currently assuming a value of '3' for impropers.")
# Permute non-central atoms
others = [atom_indices[0], atom_indices[2], atom_indices[3]]
# ((0, 1, 2), (1, 2, 0), and (2, 0, 1)) are the three paths around the trefoil
for p in [(others[i], others[j], others[k]) for (i, j, k) in [(0, 1, 2), (1, 2, 0), (2, 0, 1)]]:
# The torsion force gets added three times, since the k is divided by three
force.addTorsion(atom_indices[1], p[0], p[1], p[2],
improper_periodicity, improper_phase, improper_k/improper_idivf)
logger.info(
'{} impropers added, each applied in a six-fold trefoil'.format(
len(impropers)))
[docs]class vdWHandler(ParameterHandler):
"""Handle SMIRNOFF ``<vdW>`` tags
.. warning :: This API is experimental and subject to change.
"""
[docs] class vdWType(ParameterType):
"""A SMIRNOFF vdWForce type.
.. warning :: This API is experimental and subject to change.
"""
_VALENCE_TYPE = 'Atom' # ChemicalEnvironment valence type expected for SMARTS
_ELEMENT_NAME = 'Atom'
_SMIRNOFF_ATTRIBS = ['smirks', 'epsilon'] # Attributes expected per the SMIRNOFF spec.
_OPTIONAL_ATTRIBS = ['id', 'parent_id', 'sigma', 'rmin_half']
_REQUIRE_UNITS = {
'epsilon': unit.kilocalorie_per_mole,
'sigma': unit.angstrom,
'rmin_half': unit.angstrom
}
def __init__(self, **kwargs):
sigma = kwargs.get('sigma', None)
rmin_half = kwargs.get('rmin_half', None)
if (sigma is None) and (rmin_half is None):
raise SMIRNOFFSpecError("Either sigma or rmin_half must be specified.")
if (sigma is not None) and (rmin_half is not None):
raise SMIRNOFFSpecError(
"BOTH sigma and rmin_half cannot be specified simultaneously."
)
# TODO: Is it necessary to force everything to be sigma? We could handle having either when create_force runs
if (rmin_half is not None):
kwargs['sigma'] = 2. * rmin_half / (2.**(1. / 6.))
del kwargs['rmin_half']
super().__init__(**kwargs)
# @property
# def attrib(self):
# """Return all storable attributes as a dict.
# """
# names = ['smirks', 'sigma', 'epsilon']
# return {
# name: getattr(self, name)
# for name in names if hasattr(self, name)
# }
_TAGNAME = 'vdW' # SMIRNOFF tag name to process
_INFOTYPE = vdWType # info type to store
_OPENMMTYPE = openmm.NonbondedForce # OpenMM force class to create
# _KWARGS = ['ewaldErrorTolerance',
# 'useDispersionCorrection',
# 'usePbc'] # Kwargs to catch when create_force is called
_REQUIRE_UNITS = {'switch_width': unit.angstrom,
'cutoff': unit.angstrom}
_DEFAULT_SPEC_ATTRIBS = {
'potential': 'Lennard-Jones-12-6',
'combining_rules': 'Lorentz-Berthelot',
'scale12': 0.0,
'scale13': 0.0,
'scale14': 0.5,
'scale15': 1.0,
#'pme_tolerance': 1.e-5,
'switch_width': 1.0 * unit.angstroms,
'cutoff': 9.0 * unit.angstroms,
'method': 'cutoff',
}
_ATTRIBS_TO_TYPE = {'scale12': float,
'scale13': float,
'scale14': float,
'scale15': float
}
# TODO: Is this necessary? It's used in check_compatibility but could be hard-coded.
_SCALETOL = 1e-5
[docs] def __init__(self, **kwargs):
super().__init__(**kwargs)
self._validate_parameters()
# TODO: These properties are a fast hack and should be replaced by something better
@property
def potential(self):
"""The potential used to model van der Waals interactions"""
return self._potential
@potential.setter
def potential(self, other):
"""The potential used to model van der Waals interactions"""
valid_potentials = ['Lennard-Jones-12-6']
if other not in valid_potentials:
raise IncompatibleParameterError(f"Attempted to set vdW potential to {other}. Expected "
f"one of {valid_potentials}")
self._potential = other
@property
def combining_rules(self):
"""The combining_rules used to model van der Waals interactions"""
return self._combining_rules
@combining_rules.setter
def combining_rules(self, other):
"""The combining_rules used to model van der Waals interactions"""
valid_combining_ruless = ['Lorentz-Berthelot']
if other not in valid_combining_ruless:
raise IncompatibleParameterError(f"Attempted to set vdW combining_rules to {other}. Expected "
f"one of {valid_combining_ruless}")
self._method = other
@property
def method(self):
"""The method used to handle long-range van der Waals interactions"""
return self._method
@method.setter
def method(self, other):
"""The method used to handle long-range van der Waals interactions"""
valid_methods = ['cutoff', 'PME']
if other not in valid_methods:
raise IncompatibleParameterError(f"Attempted to set vdW method to {other}. Expected "
f"one of {valid_methods}")
self._method = other
@property
def cutoff(self):
"""The cutoff used for long-range van der Waals interactions"""
return self._cutoff
@cutoff.setter
def cutoff(self, other):
"""The cutoff used for long-range van der Waals interactions"""
unit_to_check = self._REQUIRE_UNITS['cutoff']
if not unit_to_check.unit_is_compatible(other.unit):
raise IncompatibleParameterError(
f"Attempted to set vdW cutoff to {other}, which is not compatible with "
f"expected unit {unit_to_check}")
self._cutoff = other
@property
def switch_width(self):
"""The switching width used for long-range van der Waals interactions"""
return self._switch_width
@switch_width.setter
def switch_width(self, other):
"""The switching width used for long-range van der Waals interactions"""
unit_to_check = self._REQUIRE_UNITS['switch_width']
if not unit_to_check.unit_is_compatible(other.unit):
raise IncompatibleParameterError(
f"Attempted to set vdW switch_width to {other}, which is not compatible with "
f"expected unit {unit_to_check}")
self._switch_width = other
def _validate_parameters(self):
"""
Checks internal attributes, raising an exception if they are configured in an invalid way.
"""
if self._scale12 != 0.0:
raise SMIRNOFFSpecError("Current OFF toolkit is unable to handle scale12 values other than 0.0. "
"Specified 1-2 scaling was {}".format(self._scale12))
if self._scale13 != 0.0:
raise SMIRNOFFSpecError("Current OFF toolkit is unable to handle scale13 values other than 0.0. "
"Specified 1-3 scaling was {}".format(self._scale13))
if self._scale15 != 1.0:
raise SMIRNOFFSpecError("Current OFF toolkit is unable to handle scale15 values other than 1.0. "
"Specified 1-5 scaling was {}".format(self._scale15))
supported_methods = ['cutoff', 'PME']
if self._method not in supported_methods:
raise SMIRNOFFSpecError("The Open Force Field toolkit currently only supports vdW method "
"values of {}. Received unsupported value "
"{}".format(supported_methods, self._method))
elif self._method == 'cutoff':
if self._cutoff is None:
raise SMIRNOFFSpecError("If vdW method is cutoff, a cutoff distance "
"must be provided")
elif self._method == 'PME':
if self._cutoff is None:
raise SMIRNOFFSpecError("If vdW method is PME, a cutoff distance "
"must be provided")
# if self._pme_tolerance is None:
# raise SMIRNOFFSpecError("If PME vdW method is selected, a pme_tolerance value must "
# "be specified.")
if self._potential != "Lennard-Jones-12-6":
raise SMIRNOFFSpecError("vdW potential set to {}. Only 'Lennard-Jones-12-6' is currently "
"supported".format(self._potential))
if self._combining_rules != "Lorentz-Berthelot":
raise SMIRNOFFSpecError("vdW combining_rules set to {}. Only 'Lorentz-Berthelot' is currently "
"supported".format(self._combining_rules))
# TODO: Find a better way to set defaults
# TODO: Validate these values against the supported output types (openMM force kwargs?)
# TODO: Add conditional logic to assign NonbondedMethod and check compatibility
[docs] def check_handler_compatibility(self,
handler_kwargs,
assume_missing_is_default=True):
"""
Checks if a set of kwargs used to create a ParameterHandler are compatible with this ParameterHandler. This is
called if a second handler is attempted to be initialized for the same tag. If no value is given for a field, it
will be assumed to expect the ParameterHandler class default.
Parameters
----------
handler_kwargs : dict
The kwargs that would be used to construct a ParameterHandler
assume_missing_is_default : bool
If True, will assume that parameters not specified in handler_kwargs would have been set to the default.
Therefore, an exception is raised if the ParameterHandler is incompatible with the default value for a
unspecified field.
Raises
------
IncompatibleParameterError if handler_kwargs are incompatible with existing parameters.
"""
compare_attr_to_kwargs = {
self._scale12: 'scale12',
self._scale13: 'scale13',
self._scale14: 'scale14',
self._scale15: 'scale15'
}
for attr, kwarg_key in compare_attr_to_kwargs.items():
kwarg_val = handler_kwargs.get(kwarg_key,
self._DEFAULT_SPEC_ATTRIBS[kwarg_key])
if abs(kwarg_val - attr) > self._SCALETOL:
raise IncompatibleParameterError(
"Difference between '{}' values is beyond allowed tolerance {}. "
"(handler value: {}, incompatible valie: {}".format(
kwarg_key, self._SCALETOL, attr, kwarg_val))
# TODO: Test for other possible incompatibilities here -- Probably just check for string equality for now,
# detailed check will require some openMM/MD expertise)
#self._potential: 'potential',
#self._combining_rules: 'combining_rules',
#self._switch: 'switch',
#self._cutoff: 'cutoff',
#self._method:'method'
#}
def create_force(self, system, topology, **kwargs):
self._validate_parameters()
force = openmm.NonbondedForce()
# If we're using PME, then the only possible openMM Nonbonded type is LJPME
if self._method == 'PME':
# If we're given a nonperiodic box, we always set NoCutoff. Later we'll add support for CutoffNonPeriodic
if (topology.box_vectors is None):
force.setNonbondedMethod(openmm.NonbondedForce.NoCutoff)
# if (topology.box_vectors is None):
# raise SMIRNOFFSpecError("If vdW method is PME, a periodic Topology "
# "must be provided")
else:
force.setNonbondedMethod(openmm.NonbondedForce.LJPME)
force.setCutoffDistance(9. * unit.angstrom)
force.setEwaldErrorTolerance(1.e-4)
# If method is cutoff, then we currently support openMM's PME for periodic system and NoCutoff for nonperiodic
elif self._method == 'cutoff':
# If we're given a nonperiodic box, we always set NoCutoff. Later we'll add support for CutoffNonPeriodic
if (topology.box_vectors is None):
force.setNonbondedMethod(openmm.NonbondedForce.NoCutoff)
else:
force.setNonbondedMethod(openmm.NonbondedForce.PME)
force.setUseDispersionCorrection(True)
force.setCutoffDistance(self._cutoff)
system.addForce(force)
# Iterate over all defined Lennard-Jones types, allowing later matches to override earlier ones.
atoms = self.find_matches(topology)
# Create all particles.
for particle in topology.topology_particles:
force.addParticle(0.0, 1.0, 0.0)
# Set the particle Lennard-Jones terms.
for atom_key, ljtype in atoms.items():
atom_idx = atom_key[0]
force.setParticleParameters(atom_idx, 0.0, ljtype.sigma,
ljtype.epsilon)
# Check that no atoms (n.b. not particles) are missing force parameters.
self._check_all_valence_terms_assigned(assigned_terms=atoms,
valence_terms=topology.topology_atoms)
# TODO: Can we express separate constraints for postprocessing and normal processing?
[docs] def postprocess_system(self, system, topology, **kwargs):
# Create exceptions based on bonds.
# TODO: This postprocessing must occur after the ChargeIncrementModelHandler
# QUESTION: Will we want to do this for *all* cases, or would we ever want flexibility here?
bond_particle_indices = []
for bond in topology.topology_bonds:
topology_atoms = [atom for atom in bond.atoms]
bond_particle_indices.append(
(topology_atoms[0].topology_particle_index,
topology_atoms[1].topology_particle_index))
for force in system.getForces():
# TODO: Should we just store which `Force` object we are adding to and use that instead,
# to prevent interference with other kinds of forces in the future?
# TODO: Can we generalize this to allow for `CustomNonbondedForce` implementations too?
if isinstance(force, openmm.NonbondedForce):
#nonbonded.createExceptionsFromBonds(bond_particle_indices, self.coulomb14scale, self.lj14scale)
# TODO: Don't mess with electrostatic scaling here. Have a separate electrostatics handler.
force.createExceptionsFromBonds(bond_particle_indices, 0.83333,
self._scale14)
#force.createExceptionsFromBonds(bond_particle_indices, self.coulomb14scale, self._scale14)
[docs]class ElectrostaticsHandler(ParameterHandler):
"""Handles SMIRNOFF ``<Electrostatics>`` tags.
.. warning :: This API is experimental and subject to change.
"""
_TAGNAME = 'Electrostatics'
_OPENMMTYPE = openmm.NonbondedForce
_DEPENDENCIES = [vdWHandler]
_DEFAULT_SPEC_ATTRIBS = {
'method': 'PME',
'scale12': 0.0,
'scale13': 0.0,
'scale14': 0.833333,
'scale15': 1.0,
#'pme_tolerance': 1.e-5,
#'switch_width': 8.0 * unit.angstrom, # OpenMM can't support an electrostatics switch
'switch_width': 0.0 * unit.angstrom,
'cutoff': 9.0 * unit.angstrom
}
_ATTRIBS_TO_TYPE = {'scale12': float,
'scale13': float,
'scale14': float,
'scale15': float
}
_OPTIONAL_SPEC_ATTRIBS = ['cutoff', 'switch_width']
[docs] def __init__(self, **kwargs):
super().__init__(**kwargs)
self._validate_parameters()
@property
def method(self):
"""The method used to model long-range electrostatic interactions"""
return self._method
@method.setter
def method(self, other):
"""The method used to model long-range electrostatic interactions"""
valid_methods = ['PME', 'Coulomb', 'reaction-field']
if other not in valid_methods:
raise IncompatibleParameterError(f"Attempted to set electrostatics method to {other}. Expected "
f"one of {valid_methods}")
self._method = other
@property
def cutoff(self):
"""The cutoff used for long-range van der Waals interactions"""
return self._cutoff
@cutoff.setter
def cutoff(self, other):
"""The cutoff used for long-range van der Waals interactions"""
unit_to_check = self._REQUIRE_UNITS['cutoff']
if not unit_to_check.unit_is_compatible(other.unit):
raise IncompatibleParameterError(
f"Attempted to set vdW cutoff to {other}, which is not compatible with "
f"expected unit {unit_to_check}")
self._cutoff = other
@property
def switch_width(self):
"""The switching width used for long-range electrostatics interactions"""
return self._switch_width
@switch_width.setter
def switch_width(self, other):
"""The switching width used for long-range van der Waals interactions"""
unit_to_check = self._REQUIRE_UNITS['switch_width']
if not unit_to_check.unit_is_compatible(other.unit):
raise IncompatibleParameterError(
f"Attempted to set vdW switch_width to {other}, which is not compatible with "
f"expected unit {unit_to_check}")
self._switch_width = other
def _validate_parameters(self):
"""
Checks internal attributes, raising an exception if they are configured in an invalid way.
"""
if self._scale12 != 0.0:
raise IncompatibleParameterError("Current OFF toolkit is unable to handle scale12 values other than 0.0. "
"Specified 1-2 scaling was {}".format(self._scale12))
if self._scale13 != 0.0:
raise IncompatibleParameterError("Current OFF toolkit is unable to handle scale13 values other than 0.0. "
"Specified 1-3 scaling was {}".format(self._scale13))
if self._scale15 != 1.0:
raise IncompatibleParameterError("Current OFF toolkit is unable to handle scale15 values other than 1.0. "
"Specified 1-5 scaling was {}".format(self._scale15))
supported_methods = ['PME', 'Coulomb'] # 'reaction-field'
if self._method == 'reaction-field':
raise IncompatibleParameterError('The Open Force Field toolkit does not currently support reaction-field '
'electrostatics.')
if not self._method in supported_methods:
raise IncompatibleParameterError("'method' parameter in Electrostatics tag {} is not a supported "
"option. Valid methods are {}".format(self._method, supported_methods))
if self._method == 'reaction-field' or self._method == 'PME':
if self._cutoff is None:
raise SMIRNOFFSpecError("If Electrostatics method is 'reaction-field' or 'PME', then 'cutoff' must "
"also be specified")
if self._switch_width != 0.0 * unit.angstrom:
raise IncompatibleParameterError("The current implementation of the Open Force Field toolkit can not "
"support an electrostatic switching width. Currently only `0.0 angstroms` "
"is supported (SMIRNOFF data specified {})".format(self._switch_width))
def create_force(self, system, topology, **kwargs):
existing = [system.getForce(i) for i in range(system.getNumForces())]
existing = [
f for f in existing if type(f) == openmm.NonbondedForce
]
force = existing[0]
# Among other sanity checks, this ensures that the switch value is 0.
self._validate_parameters()
# Set the nonbonded method
settings_matched = False
current_nb_method = force.getNonbondedMethod()
# First, check whether the vdWHandler set the nonbonded method to LJPME, because that means
# that electrostatics also has to be PME
if (current_nb_method == openmm.NonbondedForce.LJPME) and (self._method != 'PME'):
raise IncompatibleParameterError("In current Open Force Field toolkit implementation, if vdW "
"treatment is set to LJPME, electrostatics must also be PME "
"(electrostatics treatment currently set to {}".format(self._method))
# Then, set nonbonded methods based on method keyword
if self._method == 'PME':
# Check whether the topology is nonperiodic, in which case we always switch to NoCutoff
# (vdWHandler will have already set this to NoCutoff)
# TODO: This is an assumption right now, and a bad one. See issue #219
if topology.box_vectors is None:
assert current_nb_method == openmm.NonbondedForce.NoCutoff
settings_matched = True
# raise IncompatibleParameterError("Electrostatics handler received PME method keyword, but a nonperiodic"
# " topology. Use of PME electrostatics requires a periodic topology.")
else:
if current_nb_method == openmm.NonbondedForce.LJPME:
pass
# There's no need to check for matching cutoff/tolerance here since both are hard-coded defaults
else:
force.setNonbondedMethod(openmm.NonbondedForce.PME)
force.setCutoffDistance(9. * unit.angstrom)
force.setEwaldErrorTolerance(1.e-4)
settings_matched = True
# If vdWHandler set the nonbonded method to NoCutoff, then we don't need to change anything
elif self._method == 'Coulomb':
if topology.box_vectors is None:
# (vdWHandler will have already set this to NoCutoff)
assert current_nb_method == openmm.NonbondedForce.NoCutoff
settings_matched = True
else:
raise IncompatibleParameterError("Electrostatics method set to Coulomb, and topology is periodic. "
"In the future, this will lead to use of OpenMM's CutoffPeriodic "
"Nonbonded force method, however this is not supported in the "
"current Open Force Field toolkit.")
# If the vdWHandler set the nonbonded method to PME, then ensure that it has the same cutoff
elif self._method == 'reaction-field':
if topology.box_vectors is None:
# (vdWHandler will have already set this to NoCutoff)
assert current_nb_method == openmm.NonbondedForce.NoCutoff
settings_matched = True
else:
raise IncompatibleParameterError("Electrostatics method set to reaction-field. In the future, "
"this will lead to use of OpenMM's CutoffPeriodic or CutoffNonPeriodic"
" Nonbonded force method, however this is not supported in the "
"current Open Force Field toolkit")
if not settings_matched:
raise IncompatibleParameterError("Unable to support provided vdW method, electrostatics "
"method ({}), and topology periodicity ({}) selections. Additional "
"options for nonbonded treatment may be added in future versions "
"of the Open Force Field toolkit.".format(self._method,
topology.box_vectors is not None))
# TODO: Calculate exceptions
class ChargeIncrementModelHandler(ParameterHandler):
"""Handle SMIRNOFF ``<ChargeIncrementModel>`` tags
.. warning :: This API is experimental and subject to change.
"""
class ChargeIncrementType(ParameterType):
"""A SMIRNOFF bond charge correction type.
.. warning :: This API is experimental and subject to change.
"""
_VALENCE_TYPE = 'Bond' # ChemicalEnvironment valence type expected for SMARTS
_ELEMENT_NAME = 'ChargeIncrement'
_SMIRNOFF_ATTRIBS = ['smirks', 'chargeIncrement']
_REQUIRE_UNITS = {
'chargeIncrement': unit.elementary_charge
}
_INDEXED_ATTRIBS = ['chargeIncrement']
def __init__(self, node, parent):
super().__init__(**kwargs)
_TAGNAME = 'ChargeIncrementModel' # SMIRNOFF tag name to process
_INFOTYPE = ChargeIncrementType # info type to store
_OPENMMTYPE = openmm.NonbondedForce # OpenMM force class to create or utilize
# TODO: The structure of this is still undecided
_KWARGS = ['charge_from_molecules']
_DEFAULTS = {'number_of_conformers': 10,
'quantum_chemical_method': 'AM1',
'partial_charge_method': 'CM2'}
_ALLOWED_VALUES = {'quantum_chemical_method': ['AM1'],
'partial_charge_method': ['CM2']}
def __init__(self, **kwargs):
raise NotImplementedError("ChangeIncrementHandler is not yet implemented, pending finalization of the "
"SMIRNOFF spec")
super().__init__(**kwargs)
if number_of_conformers is None:
self._number_of_conformers = self._DEFAULTS['number_of_conformers']
elif type(number_of_conformers) is str:
self._number_of_conformers = int(number_of_conformers)
else:
self._number_of_conformers = number_of_conformers
if quantum_chemical_method is None:
self._quantum_chemical_method = self._DEFAULTS['quantum_chemical_method']
elif number_of_conformers in self._ALLOWED_VALUES['quantum_chemical_method']:
self._number_of_conformers = number_of_conformers
if partial_charge_method is None:
self._partial_charge_method = self._DEFAULTS['partial_charge_method']
elif partial_charge_method in self._ALLOWED_VALUES['partial_charge_method']:
self._partial_charge_method = partial_charge_method
def check_handler_compatibility(self, handler_kwargs, assume_missing_is_default=True):
"""
Checks if a set of kwargs used to create a ParameterHandler are compatible with this ParameterHandler. This is
called if a second handler is attempted to be initialized for the same tag. If no value is given for a field, it
will be assumed to expect the ParameterHandler class default.
Parameters
----------
handler_kwargs : dict
The kwargs that would be used to construct a ParameterHandler
assume_missing_is_default : bool
If True, will assume that parameters not specified in handler_kwargs would have been set to the default.
Therefore, an exception is raised if the ParameterHandler is incompatible with the default value for a
unspecified field.
Raises
------
IncompatibleParameterError if handler_kwargs are incompatible with existing parameters.
"""
compare_kwarg_to_attr = {
'number_of_conformers': self._number_of_conformers,
'quantum_chemical_method': self._quantum_chemical_method,
'partial_charge_method': self._partial_charge_method,
}
for kwarg_key, attr in compare_kwarg_to_attr.items():
# Skip this comparison if the kwarg isn't in handler_kwargs and we're not comparing against defaults
if not(assume_missing_is_default) and not(kwarg_key in handler_kwargs):
continue
kwarg_val = handler_kwargs.get(kwarg_key, self._DEFAULTS[kwarg_key])
if kwarg_val != attr:
raise IncompatibleParameterError(
"Incompatible '{}' values found during handler compatibility check."
"(existing handler value: {}, new existing value: {}".format(kwarg_key, attr, kwarg_val))
def assign_charge_from_molecules(self, molecule, charge_mols):
"""
Given an input molecule, checks against a list of molecules for an isomorphic match. If found, assigns
partial charges from the match to the input molecule.
Parameters
----------
molecule : an openforcefield.topology.FrozenMolecule
The molecule to have partial charges assigned if a match is found.
charge_mols : list of [openforcefield.topology.FrozenMolecule]
A list of molecules with charges already assigned.
Returns
-------
match_found : bool
Whether a match was found. If True, the input molecule will have been modified in-place.
"""
from networkx.algorithms.isomorphism import GraphMatcher
# Define the node/edge attributes that we will use to match the atoms/bonds during molecule comparison
node_match_func = lambda x, y: ((x['atomic_number'] == y['atomic_number']) and
(x['stereochemistry'] == y['stereochemistry']) and
(x['is_aromatic'] == y['is_aromatic'])
)
edge_match_func = lambda x, y: ((x['order'] == y['order']) and
(x['stereochemistry'] == y['stereochemistry']) and
(x['is_aromatic'] == y['is_aromatic'])
)
# Check each charge_mol for whether it's isomorphic to the input molecule
for charge_mol in charge_mols:
if molecule.is_isomorphic(charge_mol):
# Take the first valid atom indexing map
ref_mol_G = molecule.to_networkx()
charge_mol_G = charge_mol.to_networkX()
GM = GraphMatcher(
charge_mol_G,
ref_mol_G,
node_match=node_match_func,
edge_match=edge_match_func)
for mapping in GM.isomorphisms_iter():
topology_atom_map = mapping
break
# Set the partial charges
charge_mol_charges = charge_mol.get_partial_charges()
temp_mol_charges = charge_mol_charges.copy()
for charge_idx, ref_idx in topology_atom_map:
temp_mol_charges[ref_idx] = charge_mol_charges[charge_idx]
molecule.set_partial_charges(temp_mol_charges)
return True
# If no match was found, return False
return False
def create_force(self, system, topology, **kwargs):
from openforcefield.topology import FrozenMolecule, TopologyAtom, TopologyVirtualSite
existing = [system.getForce(i) for i in range(system.getNumForces())]
existing = [f for f in existing if type(f) == self._OPENMMTYPE]
if len(existing) == 0:
force = self._OPENMMTYPE()
system.addForce(force)
else:
force = existing[0]
for ref_mol in topology.reference_molecules:
# Make a temporary copy of ref_mol to assign charges from charge_mol
temp_mol = FrozenMolecule(ref_mol)
# First, check whether any of the reference molecules in the topology are in the charge_from_mol list
charges_from_charge_mol = False
if 'charge_from_mol' in kwargs:
charges_from_charge_mol = self.assign_charge_from_molecules(temp_mol, kwargs['charge_from_mol'])
# If the molecule wasn't assigned parameters from a manually-input charge_mol, calculate them here
if not(charges_from_charge_mol):
temp_mol.generate_conformers(num_conformers=10)
temp_mol.compute_partial_charges(quantum_chemical_method=self._quantum_chemical_method,
partial_charge_method=self._partial_charge_method)
# Assign charges to relevant atoms
for topology_molecule in topology._reference_molecule_to_topology_molecules[ref_mol]:
for topology_particle in topology_molecule.particles:
topology_particle_index = topology_particle.topology_particle_index
if type(topology_particle) is TopologyAtom:
ref_mol_particle_index = topology_particle.atom.molecule_particle_index
if type(topology_particle) is TopologyVirtualSite:
ref_mol_particle_index = topology_particle.virtual_site.molecule_particle_index
particle_charge = temp_mol._partial_charges[ref_mol_particle_index]
# Retrieve nonbonded parameters for reference atom (charge not set yet)
_, sigma, epsilon = force.getParticleParameters(topology_particle_index)
# Set the nonbonded force with the partial charge
force.setParticleParameters(topology_particle_index,
particle_charge, sigma,
epsilon)
# TODO: Move chargeModel and library residue charges to SMIRNOFF spec
def postprocess_system(self, system, topology, **kwargs):
bonds = self.find_matches(topology)
# Apply bond charge increments to all appropriate force groups
# QUESTION: Should we instead apply this to the Topology in a preprocessing step, prior to spreading out charge onto virtual sites?
for force in system.getForces():
if force.__class__.__name__ in [
'NonbondedForce'
]: # TODO: We need to apply this to all Force types that involve charges, such as (Custom)GBSA forces and CustomNonbondedForce
for (atoms, bond) in bonds.items():
# Get corresponding particle indices in Topology
particle_indices = tuple(
[atom.particle_index for atom in atoms])
# Retrieve parameters
[charge0, sigma0, epsilon0] = force.getParticleParameters(
particle_indices[0])
[charge1, sigma1, epsilon1] = force.getParticleParameters(
particle_indices[1])
# Apply bond charge increment
charge0 -= bond.increment
charge1 += bond.increment
# Update charges
force.setParticleParameters(particle_indices[0], charge0,
sigma0, epsilon0)
force.setParticleParameters(particle_indices[1], charge1,
sigma1, epsilon1)
# TODO: Calculate exceptions
class GBSAParameterHandler(ParameterHandler):
"""Handle SMIRNOFF ``<GBSAParameterHandler>`` tags
.. warning :: This API is experimental and subject to change.
"""
# TODO: Differentiate between global and per-particle parameters for each model.
# Global parameters for surface area (SA) component of model
SA_expected_parameters = {
'ACE': ['surface_area_penalty', 'solvent_radius'],
None: [],
}
# Per-particle parameters for generalized Born (GB) model
GB_expected_parameters = {
'HCT': ['radius', 'scale'],
'OBC1': ['radius', 'scale'],
'OBC2': ['radius', 'scale'],
}
class GBSAType(ParameterType):
"""A SMIRNOFF GBSA type.
.. warning :: This API is experimental and subject to change.
"""
_VALENCE_TYPE = 'Atom'
_ELEMENT_NAME = 'Atom' # TODO: This isn't actually in the spec
_SMIRNOFF_ATTRIBS = ['smirks', 'radius', 'scale']
_REQUIRE_UNITS = {'radius': unit.angstrom}
_ATTRIBS_TO_TYPE = {'scale': float}
def __init__(self, **kwargs):
super().__init__(**kwargs)
# # Store model parameters.
# gb_model = parent.attrib['gb_model']
# expected_parameters = GBSAParameterHandler.GB_expected_parameters[
# gb_model]
# provided_parameters = list()
# missing_parameters = list()
# for name in expected_parameters:
# if name in node.attrib:
# provided_parameters.append(name)
# value = _extract_quantity_from_xml_element(
# node, parent, name)
# setattr(self, name, value)
# else:
# missing_parameters.append(name)
# if len(missing_parameters) > 0:
# msg = 'GBSAForce: missing per-atom parameters for tag %s' % str(
# node)
# msg += 'model "%s" requires specification of per-atom parameters %s\n' % (
# gb_model, str(expected_parameters))
# msg += 'provided parameters : %s\n' % str(provided_parameters)
# msg += 'missing parameters: %s' % str(missing_parameters)
# raise Exception(msg)
# TODO: Finish this
_TAGNAME = 'GBSA'
_INFOTYPE = GBSAType
#_OPENMMTYPE =
def __init__(self, **kwargs):
super().__init__(**kwargs)
# TODO: Fix this
def parseElement(self):
# Initialize GB model
gb_model = element.attrib['gb_model']
valid_GB_models = GBSAParameterHandler.GB_expected_parameters.keys()
if not gb_model in valid_GB_models:
raise Exception(
'Specified GBSAForce model "%s" not one of valid models: %s' %
(gb_model, valid_GB_models))
self.gb_model = gb_model
# Initialize SA model
sa_model = element.attrib['sa_model']
valid_SA_models = GBSAParameterHandler.SA_expected_parameters.keys()
if not sa_model in valid_SA_models:
raise Exception(
'Specified GBSAForce SA_model "%s" not one of valid models: %s'
% (sa_model, valid_SA_models))
self.sa_model = sa_model
# Store parameters for GB and SA models
# TODO: Deep copy?
self.parameters = element.attrib
# TODO: Generalize this to allow forces to know when their OpenMM Force objects can be combined
def checkCompatibility(self, Handler):
"""
Check compatibility of this Handler with another Handlers.
"""
Handler = existing[0]
if (Handler.gb_model != self.gb_model):
raise ValueError(
'Found multiple GBSAForce tags with different GB model specifications'
)
if (Handler.sa_model != self.sa_model):
raise ValueError(
'Found multiple GBSAForce tags with different SA model specifications'
)
# TODO: Check other attributes (parameters of GB and SA models) automatically?
def create_force(self, system, topology, **args):
# TODO: Rework this
from openforcefield.typing.engines.smirnoff import gbsaforces
force_class = getattr(gbsaforces, self.gb_model)
force = force_class(**self.parameters)
system.addForce(force)
# Add all GBSA terms to the system.
expected_parameters = GBSAParameterHandler.GB_expected_parameters[
self.gb_model]
# Create all particles with parameters set to zero
atoms = self.getMatches(topology)
nparams = 1 + len(expected_parameters) # charge + GBSA parameters
params = [0.0 for i in range(nparams)]
for particle in topology.topology_particles():
force.addParticle(params)
# Set the GBSA parameters (keeping charges at zero for now)
for (atoms, gbsa_type) in atoms.items():
atom = atoms[0]
# Set per-particle parameters for assigned parameters
params = [atom.charge] + [
getattr(gbsa_type, name) for name in expected_parameters
]
force.setParticleParameters(atom.particle_index, params)