"""
This module contains the class for storing and creating/converting/writing
OpenMM-style ffxml files defining a force field
Author(s): Jason Swails
"""
from __future__ import absolute_import, print_function, division
from copy import copy as _copy
import math
from functools import wraps
from contextlib import closing
import datetime
from ..charmm.parameters import CharmmImproperMatchingMixin
from ..constants import DEFAULT_ENCODING
from ..formats.registry import FileFormatType
from ..modeller.residue import ResidueTemplate, PatchTemplate
from ..parameters import ParameterSet
from ..periodic_table import Element
from ..topologyobjects import NoUreyBradley
from .. import unit as u
from ..utils.io import genopen
from ..utils.six import add_metaclass, string_types, iteritems
from ..utils.six.moves import range
import warnings
from ..exceptions import ParameterWarning
from itertools import product, chain
from ..topologyobjects import (DihedralType, ImproperType, DrudeAtom)
from collections import OrderedDict
_have_lxml = False
try:
from lxml import etree
_have_lxml = True
except ImportError:
from xml.etree import ElementTree as etree
def _pretty_print_lxml(tree):
return etree.tostring(tree, encoding=DEFAULT_ENCODING, pretty_print=True).decode('utf-8')
def _pretty_print_xml_stdlib(tree):
from xml.dom import minidom
xml = etree.tostring(tree.getroot(), encoding=DEFAULT_ENCODING).decode('utf-8')
return minidom.parseString(xml).toprettyxml(indent=" ")
if _have_lxml:
pretty_print = _pretty_print_lxml
else:
pretty_print = _pretty_print_xml_stdlib
import logging
LOGGER = logging.getLogger(__name__)
[docs]def needs_lxml(func):
"""
Decorator to raise an ImportError if a function requires lxml but it is not
present
"""
@wraps(func)
def wrapper(*args, **kwargs):
if etree is None:
raise ImportError('required package lxml could not be found')
return func(*args, **kwargs)
return wrapper
[docs]@add_metaclass(FileFormatType)
class OpenMMParameterSet(ParameterSet, CharmmImproperMatchingMixin):
""" Class storing parameters from an OpenMM parameter set
Parameters
----------
filenames : str, list of str, file-like, or list of file-like; optional
Either the name of a file or a list of filenames from which parameters
should be parsed.
Notes
-----
Order is important in the list of files provided. The parameters are loaded
in the order they are provided, and any parameters that are specified in
multiple places are overwritten (that is, the *last* occurrence is the
parameter type that is used)
See Also
--------
:class:`parmed.parameters.ParameterSet`
"""
def __init__(self, *filenames):
super(OpenMMParameterSet, self).__init__()
if filenames:
raise NotImplementedError('Cannot yet read OpenMM Parameter sets')
self.unique_atom_types = False
@classmethod
def _remediate_residue_template(cls, params, residue):
"""
Modify non-compliant residue templates to conform with OpenMM requirements.
* To correctly detect waters, OpenMM ffxml water models must not contain
non-chemical bond constraints. Theses are removed when importing
foreign parameter sets (e.g., CHARMM) into OpenMM parameter sets,
and not restored on conversion from OpenMM to other formats
Parameters
----------
params : :class:`parmed.parameters.ParameterSet`
ParameterSet containing the list of parameters to be converted to a
OpenMM-compatible parameter set
residue : :class:`parmed.modeller.Residue`
The residue to remediate
Returns
-------
missing_parameters : bool
If True, the residue template is missing some parameters
"""
# Populate atomic numbers in residue template
# TODO: This can be removed if the parameter readers are guaranteed to populate this correctly
for atom in residue.atoms:
if atom.type not in params.atom_types_str:
warnings.warn('Residue {} contains atom type {} not found in parameter set and will be dropped.'.format(residue.name, atom.type))
return False
atom.atomic_number = params.atom_types_str[atom.type].atomic_number
# CHARMM Drude force field lists all lone pairs as being hydrogens??? Fix them.
types = dict((atom.name, atom.type) for atom in residue.atoms)
for lonepair in residue.lonepairs:
lp_atom = lonepair[1]
params.atom_types[types[lp_atom]].atomic_number = 0
residue[lp_atom].atomic_number = 0
# CHARMM Drude force field includes bonds to lone pairs. Delete them.
# The call to list() makes a copy of the list, so we don't modify a list
# we're iterating over.
for bond in list(residue.bonds):
if (bond.atom1.atomic_number == 0) or (bond.atom2.atomic_number == 0):
LOGGER.debug('Deleting bonds to virtual sites in residue {}'.format(residue.name))
residue.delete_bond(bond)
# Check waters
if residue.empirical_chemical_formula == 'H2O':
# Remove any H-H bonds if they are present
for bond in list(residue.bonds):
if (bond.atom1.element_name == 'H') and (bond.atom2.element_name == 'H'):
# Remove nonphysical H-H bonds
LOGGER.debug('Deleting H-H bond from water residue {}'.format(residue.name))
residue.delete_bond(bond)
else:
LOGGER.debug('keeping %s to %s %s' %(bond.atom1, bond.atom2, bond.atom2.element_name))
return True
[docs] @classmethod
def from_parameterset(cls, params, copy=False, remediate_residues=True, unique_atom_types=False):
"""
Instantiates an OpenMMParameterSet from another ParameterSet (or
subclass). The main thing this feature is responsible for is converting
lower-case atom type names into all upper-case and decorating the name
to ensure each atom type name is unique.
Warning
-------
Converting parameter sets to OpenMM can be lossy, and can modify the
original parameter set unless ``copy=True``:
* To correctly detect waters, OpenMM ffxml water models must not contain
non-chemical bond constraints. Theses are removed when importing
foreign parameter sets (e.g., CHARMM) into OpenMM parameter sets,
and not restored on conversion from OpenMM to other formats.
Parameters
----------
params : :class:`parmed.parameters.ParameterSet`
ParameterSet containing the list of parameters to be converted to as
OpenMM-compatible parameter set
copy : bool, optional, default=False
If True, the returned parameter set is a deep copy of ``params``. If
False, the returned parameter set is a shallow copy. Default False.
remediate_residues : bool, optional, default=True
If True, will remove non-chemical bonds and drop Residue definitions
that are missing parameters
unique_atom_types : bool
If True, a unique OpenMM atom type will be created for every atom of
every residue. In this case, the :class:`AtomType` objects correspond
to atom classes rather than atom types (in the OpenMM terminology).
Returns
-------
new_params : OpenMMParameterSet
OpenMMParameterSet with the same parameters as that defined in the
input parameter set
"""
new_params = cls()
if copy:
# Make a copy so we don't modify the original
params = _copy(params)
new_params.atom_types = new_params.atom_types_str = params.atom_types
new_params.atom_types_int = params.atom_types_int
new_params.atom_types_tuple = params.atom_types_tuple
new_params.bond_types = params.bond_types
new_params.angle_types = params.angle_types
new_params.urey_bradley_types = params.urey_bradley_types
new_params.dihedral_types = params.dihedral_types
new_params.improper_types = params.improper_types
new_params.improper_periodic_types = params.improper_periodic_types
new_params.rb_torsion_types = params.rb_torsion_types
new_params.cmap_types = params.cmap_types
new_params.nbfix_types = params.nbfix_types
new_params.pair_types = params.pair_types
new_params.parametersets = params.parametersets
new_params._combining_rule = params.combining_rule
new_params.default_scee = params.default_scee
new_params.default_scnb = params.default_scnb
new_params.unique_atom_types = unique_atom_types
# Copy CHARMM improper type map, if present, since this is needed for matching impropers
if hasattr(params, '_improper_key_map'):
new_params._improper_key_map = new_params._improper_key_map
if remediate_residues:
# Add only ResidueTemplate instances (no ResidueTemplateContainers)
# Maintain original residue ordering
remediated_residues = list()
for name, residue in iteritems(params.residues):
if isinstance(residue, ResidueTemplate):
# Don't discard the residue, but fix it if we need to
cls._remediate_residue_template(new_params, residue)
remediated_residues.append(residue)
for residue in remediated_residues:
new_params.residues[residue.name] = residue
for name, patch in iteritems(params.patches):
cls._remediate_residue_template(new_params, patch)
new_params.patches[patch.name] = patch
else:
# Don't remediate residues; just copy
for name, residue in iteritems(params.residues):
new_params.residues[residue.name] = residue
for name, patch in iteritems(params.patches):
new_params.patches[patch.name] = patch
# Only add unique patches
unique_patches = OrderedDict()
discarded_patches = []
for name, patch in iteritems(params.patches):
if isinstance(patch, PatchTemplate):
templhash = OpenMMParameterSet._templhasher(patch)
if templhash not in unique_patches:
new_params.patches[name] = patch
unique_patches[templhash] = patch
else:
patch_collision = unique_patches[templhash]
warnings.warn('Patch {} discarded because OpenMM considers it identical to {}'.format(patch, patch_collision))
discarded_patches.append(patch)
if (len(discarded_patches) > 0):
warnings.warn('{} patches discarded, {} retained'.format(len(discarded_patches), len(new_params.patches)))
for patch in discarded_patches:
del new_params.patches[patch.name]
return new_params
def _get_mm_atom_type(self, atom, residue, drude=False):
"""Get the OpenMM atom type for an atom.
Parameters
----------
atom : :class:`Atom`
the atom for which to get the type
residue : :class:`ResidueTemplate` or :class:`PatchTemplate`
the residue the atom belongs to
drude : bool
if True, get the atom type for the Drude particle attached to the
atom rather than the atom itself
"""
if self.unique_atom_types:
if drude:
return 'Drude-%s-%s' % (residue.name, atom.name)
return '%s-%s' % (residue.name, atom.name)
if drude:
return atom.drude_type
return atom.type
[docs] @needs_lxml
def write(self, dest, provenance=None, write_unused=True, separate_ljforce=False,
improper_dihedrals_ordering='default', charmm_imp=False,
skip_duplicates=True):
""" Write the parameter set to an XML file for use with OpenMM
Parameters
----------
dest : str or file-like
The name of the file or the file-like object (with a ``write``
attribute) to which the XML file will be written
provenance : dict, optional
If present, the XML file will be tagged with the available fields.
Keys of the dictionary become XML etree.Element tags, the values of the
dictionary must be instances of any of:
- str / unicode (Py2) or str (Py3) - one XML element with this
content is written
- list - one XML element per each item of the list is written, all
these XML elements use the same tag (key in provenance dict)
- dict - one of the keys of this dict must be the same as the key of
of the provenance dict under which this dict is nested. The value
under this key becomes the content of the XML element. Remaining keys
and their values are used to construct attributes of the XML element.
Note that OrderedDict's should be used to ensure appropriate order
of the XML elements and their attributes.
Default is no provenance.
Example (unordered):
provenance = {'Reference' : ['Nature', 'Cell'],
'Source' : {'Source': 'leaprc.ff14SB', sourcePackage :
'AmberTools', sourcePackageVersion : '15'},
'User' : 'Mark'}
write_unused : bool
If False: a) residue templates using unavailable atom types will not
be written, b) atom types that are not used in any of the residue
templates remaining and parameters including those atom types will
not be written. A ParameterWarning is issued if any such residues are
found in a).
separate_ljforce : bool
If True will use a separate LennardJonesForce to create a
CostumNonbondedForce to compute L-J interactions. It will set sigma
to 1 and epsilon to 0 in the NonbondedForce so that the
NonbondedForce only calculates the electrostatic contribution. It
should be set to True when converting a CHARMM force field file that
doesn't have pair-specific L-J modifications (NBFIX in CHARMM) so
that the ffxml conversion is compatible with the main charmm36.xml file.
Note:
----
When pair-specific L-J modifications are present (NBFIX in CHARMM), this
behavior is always present and this flag is ignored.
improper_dihedrals_ordering : str
The ordering to use when assigning improper torsions in OpenMM. Default is 'default',
other option is 'amber'
charmm_imp: bool
If True, will check for existence of IMPR in each residue and patch template,
and write out the explicit improper definition without wildcards in the ffxml file.
skip_duplicates : bool
If True: residues which appear identical to an existing residue will
not be written. This is usually the best choice. The most common
reason for setting skip_duplicates to false is if you have different
parametrizations for stereoisomers. Note that since OpenMM's residue
hashing is not aware of chirality, if you wish to use the results in
simulations you will need to explicitly provide the template names
for affected residues.
Notes
-----
The generated XML file will have the XML tag ``DateGenerated`` added to
the provenance information set to the current date. Therefore, you
should not provide this information in ``provenance`` (it will be
removed if it is provided).
"""
if not write_unused:
skip_residues = self._find_unused_residues()
skip_types = self._find_unused_types(skip_residues)
if skip_residues:
warnings.warn('Some residue templates using unavailable AtomTypes '
'were found. They will not be written to the ffxml '
'as write_unused is set to False', ParameterWarning)
else:
skip_residues = set()
skip_types = set()
if self.atom_types:
try:
self.typeify_templates()
except KeyError:
warnings.warn('Some residue templates are using unavailable '
'AtomTypes', ParameterWarning)
[valid_residues_for_patch, valid_patches_for_residue] = self._determine_valid_patch_combinations(skip_residues)
LOGGER.debug('Valid patch combinations:')
for patch_name in self.patches:
LOGGER.debug('%8s : %s', patch_name, valid_residues_for_patch[patch_name])
if charmm_imp:
self._find_explicit_impropers()
self._compress_impropers()
root = etree.Element('ForceField')
self._write_omm_provenance(root, provenance)
self._write_omm_atom_types(root, skip_types, skip_residues)
self._write_omm_residues(root, skip_residues, skip_duplicates,
valid_patches_for_residue=valid_patches_for_residue)
self._write_omm_patches(root, valid_residues_for_patch)
self._write_omm_bonds(root, skip_types)
self._write_omm_angles(root, skip_types)
self._write_omm_urey_bradley(root, skip_types)
self._write_omm_dihedrals(root, skip_types, improper_dihedrals_ordering)
self._write_omm_impropers(root, skip_types)
#self._write_omm_rb_torsions(root, skip_types)
self._write_omm_cmaps(root, skip_types)
self._write_omm_scripts(root, skip_types)
self._write_omm_nonbonded(root, skip_types, separate_ljforce)
self._write_omm_LennardJonesForce(root, skip_types, separate_ljforce)
self._write_omm_DrudeForce(root, skip_types)
tree = etree.ElementTree(root)
xml = pretty_print(tree)
if isinstance(dest, string_types):
with closing(genopen(dest, 'w')) as f:
f.write(xml)
else:
dest.write(xml)
def _find_explicit_impropers(self):
"""
For every residue, find any explicitly-specified (e.g. CHARMM) improper torsions and identify all wild-card improper parameters that match.
Expand all of these out into explicit impropers.
This is necessary for OpenMM to correctly handle impropers for these residues.
.. todo ::
* Do we need to do this for patches as well?
"""
# Regenerate improper key map
self._improper_key_map = OrderedDict()
for key in self.improper_types.keys():
self._improper_key_map[tuple(sorted(key))] = key
improper_harmonic = OrderedDict() # improper_harmonic[key] is the harmonic improper parameter for unique key `key`
improper_periodic = OrderedDict() # improper_harmonic[key] is the periodic improper parameter for unique key `key`
C_types = [t for t in self.atom_types if self.atom_types[t].atomic_number == 6]
N_types = [t for t in self.atom_types if self.atom_types[t].atomic_number == 7]
def get_types(residue, atomname):
"""Return list of atom type(s) that match the given atom name.
"""
a_names = [a.name for a in residue.atoms]
a_types = [a.type for a in residue.atoms]
if atomname == '-C':
return C_types
elif atomname == '+N':
return N_types
elif atomname[0] in ['-', '+']:
raise ValueError('Unknown atom name %s' % atomname)
else:
return [ a_types[a_names.index(atomname)] ]
# Iterate over all residues
for name, residue in chain(iteritems(self.residues), iteritems(self.patches)):
for impr in residue._impr:
# Get the list of types involved in this improper
try:
types = [ get_types(residue, atomname) for atomname in impr ]
except ValueError:
continue
improper_found = False
for key in product(*types):
# Search for an improper that matches these types
improper = self.match_improper_type(*key)
if improper is None:
continue
# Add this to our types
if isinstance(improper, ImproperType):
improper_harmonic[key] = improper
improper_found = True
elif isinstance(improper, DihedralType):
improper_periodic[key] = improper
improper_found = True
else:
raise Exception('Something went wrong with improper type for {} returning an unexpected object {}'.format(key, improper))
# Warn if no improper was found
if not improper_found:
raise Exception('No improper found for improper {} in residue {} (types were {})'.format(impr, name, types))
# Update our impropers
self.improper_periodic_types = improper_periodic
self.improper_types = improper_harmonic
def _compress_impropers(self):
"""
OpenMM's ForceField cannot handle impropers that match the same four atoms
in more than one order, so Peter Eastman wants us to compress duplicates
and increment the spring constant accordingly.
"""
if not self.improper_types: return
unique_keys = OrderedDict() # unique_keys[key] is the key to retrieve the improper from improper_types
improper_types = OrderedDict() # replacement for self.improper_types with compressed impropers
for atoms, improper in iteritems(self.improper_types):
# Compute a unique key
unique_key = tuple(sorted(atoms))
if unique_key in unique_keys:
# Accumulate spring constant, discarding this contribution
# TODO: Do we need to check if `psi_eq` is the same?
atoms2 = unique_keys[unique_key]
improper_types[atoms2].psi_k += improper.psi_k
warnings.warn('Compressing improper {} because it contains same atoms as {}'.format(improper, improper_types[atoms2]))
else:
# Store this improper
unique_keys[unique_key] = atoms
improper_types[atoms] = improper
self.improper_types = improper_types
def _find_unused_residues(self):
skip_residues = set()
for name, residue in iteritems(self.residues):
if any((atom.type not in self.atom_types for atom in residue.atoms)):
skip_residues.add(name)
return skip_residues
def _find_unused_types(self, skip_residues):
keep_types = set()
for name, residue in iteritems(self.residues):
if name not in skip_residues:
for atom in residue.atoms:
keep_types.add(atom.type)
return {typ for typ in self.atom_types if typ not in keep_types}
@staticmethod
def _templhasher(residue):
"""
Create a unique hash for each residue and patch template using only properties rendered to OpenMM ffxml.
"""
hash_info = tuple()
# Sort tuples of atom properties by atom name
if len(residue.atoms) > 0:
hash_info += tuple(sorted( [(atom.type, str(atom.charge)) for atom in residue.atoms] ))
# Sort list of deleted atoms by atom name
if hasattr(residue, 'delete_atoms') and len(residue.delete_atoms) > 0:
hash_info += tuple(sorted([atom_name for atom_name in residue.delete_atoms]))
# Sort list of bonds by first bond name
if len(residue.bonds) > 0:
hash_info += tuple(sorted( [(bond.atom1.name, bond.atom2.name) if (bond.atom1.name < bond.atom2.name) else (bond.atom2.name, bond.atom1.name) for bond in residue.bonds] ))
# Add head and tail
if residue.head:
hash_info += (residue.head.name,)
if residue.tail:
hash_info += (residue.tail.name,)
# TODO: Is there any other data that is rendered to ffxml files we should include?
return hash(hash_info)
@needs_lxml
def _write_omm_provenance(self, root, provenance):
info = etree.SubElement(root, 'Info')
date_generated = etree.SubElement(info, "DateGenerated")
date_generated.text = '%02d-%02d-%02d' % datetime.datetime.now().timetuple()[:3]
provenance = provenance or OrderedDict()
for tag, content in iteritems(provenance):
if tag == 'DateGenerated': continue
if not isinstance(content, list):
content = [content]
for sub_content in content:
if isinstance(sub_content, string_types):
item = etree.Element(tag)
item.text = sub_content
info.append(item)
elif isinstance(sub_content, OrderedDict) or isinstance(sub_content, dict):
if tag not in sub_content:
raise KeyError('Content of an attribute-containing element '
'specified incorrectly.')
attributes = [key for key in sub_content if key != tag]
element_content = sub_content[tag]
attributes = { k : str(v) for (k,v) in sub_content.items() }
item = etree.SubElement(info, tag, **attributes)
item.text = str(element_content)
else:
raise TypeError('Incorrect type of the %s element content' % tag)
@needs_lxml
def _write_omm_atom_types(self, xml_root, skip_types, skip_residues):
if not self.atom_types: return
xml_section = etree.SubElement(xml_root, "AtomTypes")
if self.unique_atom_types:
for residue in list(self.residues.values())+list(self.patches.values()):
if residue.name in skip_residues: continue
for atom in residue.atoms:
atom_type = self.atom_types[atom.type]
properties = { 'name' : self._get_mm_atom_type(atom, residue), 'class' : atom.type, 'mass' : str(atom_type.mass) }
if atom_type.atomic_number != 0:
properties['element'] = str(Element[atom_type.atomic_number])
etree.SubElement(xml_section, 'Type', **properties)
if isinstance(atom, DrudeAtom):
properties = { 'name' : self._get_mm_atom_type(atom, residue, True), 'class' : atom.drude_type, 'mass' : '0.0' }
etree.SubElement(xml_section, 'Type', **properties)
else:
for name, atom_type in iteritems(self.atom_types):
if name in skip_types: continue
assert atom_type.atomic_number >= 0, 'Atomic number not set!'
properties = { 'name' : name, 'class' : name, 'mass' : str(atom_type.mass) }
if atom_type.atomic_number == 0:
etree.SubElement(xml_section, 'Type', **properties)
else:
element = Element[atom_type.atomic_number]
etree.SubElement(xml_section, 'Type', element=str(element), **properties)
def _get_lonepair_parameters(self, lonepair):
(lptype, a1, a2, a3, a4, r, theta, phi) = lonepair
if lptype == 'relative':
xweights = [-1.0, 0.0, 1.0]
elif lptype == 'bisector':
xweights = [-1.0, 0.5, 0.5]
else:
raise ValueError('Unknown lonepair type: '+lptype)
r /= 10.0 # convert to nanometers
theta *= math.pi / 180.0 # convert to radians
phi = (180 - phi) * math.pi / 180.0 # convert to radians
p = [r*math.cos(theta), r*math.sin(theta)*math.cos(phi), r*math.sin(theta)*math.sin(phi)]
p = [x if abs(x) > 1e-10 else 0 for x in p] # Avoid tiny numbers caused by roundoff error
return dict(type="localCoords",
siteName=a1, atomName1=a2, atomName2=a3, atomName3=a4,
wo1="1", wo2="0", wo3="0",
wx1=str(xweights[0]), wx2=str(xweights[1]), wx3=str(xweights[2]),
wy1="0", wy2="-1", wy3="1",
p1=str(p[0]), p2=str(p[1]), p3=str(p[2]))
@needs_lxml
def _write_omm_residues(self, xml_root, skip_residues, skip_duplicates, valid_patches_for_residue=None):
if not self.residues: return
if valid_patches_for_residue is None:
valid_patches_for_residue = OrderedDict()
written_residues = OrderedDict()
xml_section = etree.SubElement(xml_root, 'Residues')
for name, residue in iteritems(self.residues):
if name in skip_residues: continue
templhash = OpenMMParameterSet._templhasher(residue)
if templhash in written_residues:
residue_collision = written_residues[templhash]
if skip_duplicates:
warnings.warn('Skipping writing of residue {} because OpenMM considers it identical to {}'.format(residue, residue_collision))
continue
else:
warnings.warn('Residue {} will be considered by OpenMM to be identical to {}.'.format(residue, residue_collision))
written_residues[templhash] = residue
# Write residue
if residue.override_level == 0:
xml_residue = etree.SubElement(xml_section, 'Residue', name=residue.name)
else:
xml_residue = etree.SubElement(xml_section, 'Residue', name=residue.name, override=str(residue.override_level))
# Write residue contents
for atom in residue.atoms:
if isinstance(atom, DrudeAtom):
etree.SubElement(xml_residue, 'Atom', name=atom.name, type=self._get_mm_atom_type(atom, residue), charge=str(atom.charge-atom.drude_charge))
etree.SubElement(xml_residue, 'Atom', name='D'+atom.name, type=self._get_mm_atom_type(atom, residue, True), charge=str(atom.drude_charge))
else:
etree.SubElement(xml_residue, 'Atom', name=atom.name, type=self._get_mm_atom_type(atom, residue), charge=str(atom.charge))
for bond in residue.bonds:
etree.SubElement(xml_residue, 'Bond', atomName1=bond.atom1.name, atomName2=bond.atom2.name)
for lonepair in residue.lonepairs:
etree.SubElement(xml_residue, 'VirtualSite', self._get_lonepair_parameters(lonepair))
for atom in residue.connections:
etree.SubElement(xml_residue, 'ExternalBond', atomName=atom.name)
if residue.head is not None:
etree.SubElement(xml_residue, 'ExternalBond', atomName=residue.head.name)
if residue.tail is not None and residue.tail is not residue.head:
etree.SubElement(xml_residue, 'ExternalBond', atomName=residue.tail.name)
if residue.name in valid_patches_for_residue:
for patch_name in valid_patches_for_residue[residue.name]:
etree.SubElement(xml_residue, 'AllowPatch', name=patch_name)
def _determine_valid_patch_combinations(self, skip_residues):
"""
Determine valid (permissible) combinations of patches with residues that
lead to integral net charges.
Parameters
----------
skip_residues : set of ResidueTemplate
List of residues to skip
Returns
-------
valid_residues_for_patch : dict
valid_residues_for_patch[patch] is a list of residues compatible with that patch
valid_patches_for_residue : dict
valid_patches_for_residue[residue] is a list of patches compatible with that residue
"""
# Attempt to patch every residue, recording only valid combinations.
valid_residues_for_patch = OrderedDict()
for patch in self.patches.values():
valid_residues_for_patch[patch.name] = list()
valid_patches_for_residue = OrderedDict()
for residue in self.residues.values():
valid_patches_for_residue[residue.name] = list()
# Create list of residues to check compatibility against
residues = [ residue for residue in self.residues.values() if (residue not in skip_residues) ]
# Check patch compatibilities
for patch in self.patches.values():
residue_compatibilities = [ residue.patch_is_compatible(patch) for residue in residues ]
for (residue, is_compatible) in zip(residues, residue_compatibilities):
if is_compatible:
valid_residues_for_patch[patch.name].append(residue.name)
valid_patches_for_residue[residue.name].append(patch.name)
return [valid_residues_for_patch, valid_patches_for_residue]
@needs_lxml
def _write_omm_patches(self, xml_root, valid_residues_for_patch, write_apply_to_residue=False):
"""
Write patch definitions for OpenMM ForceField
Parameters
----------
xml_root : lxml.etree.Element
The XML Element write the <Patches> section to.
valid_residues_for_patch : dict of str : str
valid_residues_for_patch[patch_name] lists the residue names valid for this patch
write_apply_to_residue : bool, optional, default=False
If True, will write <ApplyToResidue> tags.
"""
if not self.patches: return
written_patches = OrderedDict()
xml_patches = etree.SubElement(xml_root, 'Patches')
for name, patch in iteritems(self.patches):
# Require that at least one valid patch combination exists for this patch
if (name not in valid_residues_for_patch) or (len(valid_residues_for_patch[name])==0):
continue
templhash = OpenMMParameterSet._templhasher(patch) # TODO: this may be redundant now
if templhash in written_patches:
patch_collision = written_patches[templhash]
warnings.warn('Skipping writing of patch {} because OpenMM considers it identical to {}'.format(patch, patch_collision))
continue
written_patches[templhash] = patch
if patch.override_level == 0:
patch_xml = etree.SubElement(xml_patches, 'Patch', name=patch.name)
else:
patch_xml = etree.SubElement(xml_patches, 'Patch', name=patch.name, override=str(patch.override_level))
# To generate the patch definition, we need to apply it to a residue and see exactly what
# changes. We might get different definitions depending on which residue we pick, so try
# all possible residues to take the most common result.
versions = {}
for residue_name in valid_residues_for_patch[name]:
try:
residue = self.residues[residue_name]
except KeyError as e:
msg = 'Compatible residue not found in self.residues\n'
msg += ' patch name: %s\n' % name
msg += ' valid patch combinations: %s\n' % valid_residues_for_patch[name]
msg += ' residue name: %s\n' % residue_name
msg += str(e)
raise(msg)
patched_residue = residue.apply_patch(patch)
instructions = []
for atom in patch.atoms:
if atom.name not in residue:
command = 'AddAtom'
else:
command = 'ChangeAtom'
if isinstance(atom, DrudeAtom):
instructions.append((command, dict(name=atom.name, type=self._get_mm_atom_type(atom, patch), charge=str(atom.charge-atom.drude_charge))))
instructions.append((command, dict(name='D'+atom.name, type=self._get_mm_atom_type(atom, patch, True), charge=str(atom.drude_charge))))
else:
instructions.append((command, dict(name=atom.name, type=self._get_mm_atom_type(atom, patch), charge=str(atom.charge))))
for atom_name in patch.delete_atoms:
instructions.append(('RemoveAtom', dict(name=atom_name)))
for bond in patch.bonds:
instructions.append(('RemoveBond', dict(atomName1=bond.atom1.name, atomName2=bond.atom2.name)))
for bond in patched_residue.bonds:
if (bond.atom1.name not in residue) or (bond.atom2.name not in residue):
if (bond.atom1.atomic_number != 0) or (bond.atom2.atomic_number != 0): # CHARMM adds bonds to lone pairs, which we need to omit.
instructions.append(('AddBond', dict(atomName1=bond.atom1.name, atomName2=bond.atom2.name)))
for bond in residue.bonds:
if (bond.atom1.name not in patched_residue) or (bond.atom2.name not in patched_residue):
instructions.append(('RemoveBond', dict(atomName1=bond.atom1.name, atomName2=bond.atom2.name)))
if (residue.head is not None) and (patched_residue.head is None):
instructions.append(('RemoveExternalBond', dict(atomName=residue.head.name)))
if (residue.tail is not None) and (patched_residue.tail is None):
instructions.append(('RemoveExternalBond', dict(atomName=residue.tail.name)))
if (residue.head is None) and (patched_residue.head is not None):
instructions.append(('AddExternalBond', dict(atomName=patched_residue.head.name)))
if (residue.tail is None) and (patched_residue.tail is not None):
instructions.append(('AddExternalBond', dict(atomName=patched_residue.tail.name)))
for lonepair in patch.lonepairs:
instructions.append(('VirtualSite', self._get_lonepair_parameters(lonepair)))
if write_apply_to_residue:
for residue_name in valid_residues_for_patch[patch.name]:
instructions.append(('ApplyToResidue', dict(name=residue_name)))
# Convert to hashable types
instructions = tuple((i[0], tuple(item for item in i[1].items())) for i in instructions)
if instructions in versions:
versions[instructions] += 1
else:
versions[instructions] = 1
# Write the consensus definition.
max_count = max(versions.values())
instructions = [key for key, value in versions.items() if value == max_count][0]
for command, attrib in instructions:
etree.SubElement(patch_xml, command, dict(attrib))
@needs_lxml
def _write_omm_bonds(self, xml_root, skip_types):
if not self.bond_types: return
xml_force = etree.SubElement(xml_root, 'HarmonicBondForce')
bonds_done = set()
lconv = u.angstroms.conversion_factor_to(u.nanometers)
kconv = u.kilocalorie.conversion_factor_to(u.kilojoule) / lconv**2 * 2
for (a1, a2), bond in iteritems(self.bond_types):
if any((a in skip_types for a in (a1, a2))): continue
if (a1, a2) in bonds_done: continue
bonds_done.add((a1, a2))
bonds_done.add((a2, a1))
etree.SubElement(xml_force, 'Bond', class1=a1, class2=a2, length=str(bond.req*lconv), k=str(bond.k*kconv))
@needs_lxml
def _write_omm_angles(self, xml_root, skip_types):
if not self.angle_types: return
xml_force = etree.SubElement(xml_root, 'HarmonicAngleForce')
angles_done = set()
tconv = u.degree.conversion_factor_to(u.radians)
kconv = u.kilocalorie.conversion_factor_to(u.kilojoule) * 2
for (a1, a2, a3), angle in iteritems(self.angle_types):
if any((a in skip_types for a in (a1, a2, a3))): continue
if (a1, a2, a3) in angles_done: continue
angles_done.add((a1, a2, a3))
angles_done.add((a3, a2, a1))
etree.SubElement(xml_force, 'Angle', class1=a1, class2=a2, class3=a3, angle=str(angle.theteq*tconv), k=str(angle.k*kconv))
@needs_lxml
def _write_omm_dihedrals(self, xml_root, skip_types, improper_dihedrals_ordering):
if not self.dihedral_types and not self.improper_periodic_types: return
# In ParameterSet, dihedral_types is *always* of type DihedralTypeList.
# The from_structure method ensures that, even if the containing
# Structure has separate dihedral entries for each torsion
if improper_dihedrals_ordering == 'default':
xml_force = etree.SubElement(xml_root, 'PeriodicTorsionForce')
else:
xml_force = etree.SubElement(xml_root, 'PeriodicTorsionForce', ordering=improper_dihedrals_ordering)
diheds_done = set()
pconv = u.degree.conversion_factor_to(u.radians)
kconv = u.kilocalorie.conversion_factor_to(u.kilojoule)
def nowild(name):
return name if name != 'X' else ''
for (a1, a2, a3, a4), dihed in iteritems(self.dihedral_types):
if any((a in skip_types for a in (a1, a2, a3, a4))): continue
if (a1, a2, a3, a4) in diheds_done: continue
diheds_done.add((a1, a2, a3, a4))
diheds_done.add((a4, a3, a2, a1))
terms = OrderedDict()
for i, term in enumerate(dihed):
i += 1
terms['periodicity%d' % i] = str(term.per)
terms['phase%d' % i] = str(term.phase*pconv)
terms['k%d' % i] = str(term.phi_k*kconv)
etree.SubElement(xml_force, 'Proper', class1=nowild(a1), class2=a2, class3=a3, class4=nowild(a4), **terms)
# Now do the periodic impropers. OpenMM expects the central atom to be
# listed first. ParameterSet goes out of its way to list it third
# (consistent with Amber) except in instances where order is random (as
# in CHARMM parameter files). But CHARMM parameter files don't have
# periodic impropers, so we don't have to worry about that here.
for (a2, a3, a1, a4), improp in iteritems(self.improper_periodic_types):
if any((a in skip_types for a in (a1, a2, a3, a4))): continue
# Try to make the wild-cards in the middle
if a4 == 'X':
if a2 != 'X':
a2, a4 = a4, a2
elif a3 != 'X':
a3, a4 = a4, a3
if a2 != 'X' and a3 == 'X':
# Single wild-card entries put the wild-card in position 2
a2, a3 = a3, a2
etree.SubElement(xml_force, 'Improper', class1=a1, class2=nowild(a2), class3=nowild(a3), class4=nowild(a4),
periodicity1=str(improp.per), phase1=str(improp.phase*pconv), k1=str(improp.phi_k*kconv))
@needs_lxml
def _write_omm_impropers(self, xml_root, skip_types):
if not self.improper_types: return
xml_force = etree.SubElement(xml_root, 'CustomTorsionForce', energy="k*(theta-theta0)^2")
etree.SubElement(xml_force, 'PerTorsionParameter', name="k")
etree.SubElement(xml_force, 'PerTorsionParameter', name="theta0")
kconv = u.kilocalorie.conversion_factor_to(u.kilojoule)
tconv = u.degree.conversion_factor_to(u.radian)
def nowild(name):
return name if name != 'X' else ''
for (a1, a2, a3, a4), improp in iteritems(self.improper_types):
if any((a in skip_types for a in (a1, a2, a3, a4))): continue
etree.SubElement(xml_force, 'Improper', class1=nowild(a1), class2=nowild(a2), class3=nowild(a3), class4=nowild(a4),
k=str(improp.psi_k*kconv), theta0=str(improp.psi_eq*tconv))
@needs_lxml
def _write_omm_urey_bradley(self, xml_root, skip_types):
if not self.urey_bradley_types: return None
xml_root.append( etree.Comment("Urey-Bradley terms") )
xml_force = etree.SubElement(xml_root, 'AmoebaUreyBradleyForce')
length_conv = u.angstroms.conversion_factor_to(u.nanometers)
_ambfrc = u.kilocalorie_per_mole/u.angstrom**2
_ommfrc = u.kilojoule_per_mole/u.nanometer**2
frc_conv = _ambfrc.conversion_factor_to(_ommfrc)
ureys_done = set()
for (a1, a2, a3), urey in iteritems(self.urey_bradley_types):
if any((a in skip_types for a in (a1, a2, a3))): continue
if (a1, a2, a3) in ureys_done: continue
if urey == NoUreyBradley: continue
etree.SubElement(xml_force, 'UreyBradley', class1=a1, class2=a2, class3=a3, d=str(urey.req*length_conv), k=str(urey.k*frc_conv))
@needs_lxml
def _write_omm_cmaps(self, xml_root, skip_types):
if not self.cmap_types: return
xml_force = etree.SubElement(xml_root, 'CMAPTorsionForce')
maps = OrderedDict()
counter = 0
econv = u.kilocalorie.conversion_factor_to(u.kilojoule)
for _, cmap in iteritems(self.cmap_types):
if id(cmap) in maps: continue
maps[id(cmap)] = counter
counter += 1
xml_map = etree.SubElement(xml_force, 'Map')
grid = cmap.grid.switch_range().T
map_string = ''
for i in range(cmap.resolution):
base = i * cmap.resolution
for j in range(cmap.resolution):
map_string += ' %s' % (grid[base+j]*econv)
map_string += '\n'
xml_map.text = map_string
used_torsions = set()
for (a1, a2, a3, a4, _, _, _, a5), cmap in iteritems(self.cmap_types):
if any((a in skip_types for a in (a1, a2, a3, a4, a5))): continue
if (a1, a2, a3, a4, a5) in used_torsions: continue
used_torsions.add((a1, a2, a3, a4, a5))
used_torsions.add((a5, a4, a3, a2, a1))
etree.SubElement(xml_force, 'Torsion', map=str(maps[id(cmap)]),
class1=a1, class2=a2, class3=a3, class4=a4, class5=a5)
@needs_lxml
def _write_omm_nonbonded(self, xml_root, skip_types, separate_ljforce):
if not self.atom_types: return
# Compute conversion factors for writing in natrual OpenMM units.
length_conv = u.angstrom.conversion_factor_to(u.nanometer)
ene_conv = u.kilocalories.conversion_factor_to(u.kilojoules)
# Get the 1-4 scaling factors from the torsion list
scee, scnb = set(), set()
for key in self.dihedral_types:
dt = self.dihedral_types[key]
for t in dt:
if t.scee: scee.add(t.scee)
if t.scnb: scnb.add(t.scnb)
if len(scee) > 1:
raise NotImplementedError('Cannot currently handle mixed 1-4 '
'scaling: Elec. Scaling factors %s detected' %
(', '.join([str(x) for x in scee])))
if len(scnb) > 1:
raise NotImplementedError('Cannot currently handle mixed 1-4 '
'scaling: L-J Scaling factors %s detected' %
(', '.join([str(x) for x in scnb])))
if len(scee) > 0:
coulomb14scale = 1.0 / scee.pop()
else:
coulomb14scale = 1.0 / self.default_scee
if len(scnb) > 0:
lj14scale = 1.0 / scnb.pop()
else:
lj14scale = 1.0 / self.default_scnb
# Write NonbondedForce records.
xml_force = etree.SubElement(xml_root, 'NonbondedForce', coulomb14scale=str(coulomb14scale), lj14scale=str(lj14scale))
etree.SubElement(xml_force, 'UseAttributeFromResidue', name="charge")
for name, atom_type in iteritems(self.atom_types):
if name in skip_types: continue
if (atom_type.rmin is not None) and (atom_type.epsilon is not None):
sigma = atom_type.sigma * length_conv # in md_unit_system
epsilon = atom_type.epsilon * ene_conv # in md_unit_system
else:
# Dummy atom
sigma = 1.0
epsilon = 0.0
if self.nbfix_types or separate_ljforce:
# turn off L-J. Will use LennardJonesForce to use CostumNonbondedForce to compute L-J interactions
sigma = 1.0
epsilon = 0.0
else:
# NonbondedForce cannot handle distinct 14 parameters
# We need to use a separate LennardJonesForce instead
# TODO: Can we autodetect this and switch on separate_ljforce earlier?
if (atom_type.rmin_14 != atom_type.rmin) or (atom_type.epsilon_14 != atom_type.epsilon):
raise NotImplementedError('OpenMM <NonbondedForce> cannot handle '
'distinct 1-4 sigma and epsilon parameters; '
'use separate_ljforce=True instead')
# Ensure we don't have sigma = 0
if (sigma == 0.0):
if (epsilon == 0.0):
sigma = 1.0 # reset sigma = 1
else:
raise ValueError("For atom type '%s', sigma = 0 but "
"epsilon != 0." % name)
attributes = { 'class' : name, 'sigma' : str(sigma), 'epsilon' : str(abs(epsilon)) }
etree.SubElement(xml_force, 'Atom', **attributes)
@needs_lxml
def _write_omm_LennardJonesForce(self, xml_root, skip_types, separate_ljforce):
if not self.nbfix_types and not separate_ljforce: return
# Convert Conversion factors for writing in natural OpenMM units
length_conv = u.angstrom.conversion_factor_to(u.nanometer)
ene_conv = u.kilocalories.conversion_factor_to(u.kilojoules)
scnb = set()
for key in self.dihedral_types:
dt = self.dihedral_types[key]
for t in dt:
if t.scnb: scnb.add(t.scnb)
if len(scnb) > 1:
raise NotImplementedError('Cannot currently handle mixed 1-4 '
'scaling: L-J Scaling factors %s detected' %
(', '.join([str(x) for x in scnb])))
if len(scnb) > 0:
lj14scale = 1.0 / scnb.pop()
else:
lj14scale = 1.0 / self.default_scnb
# write L-J records
xml_force = etree.SubElement(xml_root, 'LennardJonesForce', lj14scale=str(lj14scale))
for name, atom_type in iteritems(self.atom_types):
if name in skip_types: continue
if (atom_type.rmin is not None) and (atom_type.epsilon is not None):
sigma = atom_type.sigma * length_conv # in md_unit_system
epsilon = atom_type.epsilon * ene_conv # in md_unit_system
else:
# Dummy atom
sigma = 1.0
epsilon = 0.0
# Ensure we don't have sigma = 0
if (sigma == 0.0):
if (epsilon == 0.0):
sigma = 1.0 # reset sigma = 1
else:
raise ValueError("For atom type '%s', sigma = 0 but "
"epsilon != 0." % name)
# Handle special values used for 14 interactions
if (atom_type.rmin_14 != atom_type.rmin) or (atom_type.epsilon_14 != atom_type.epsilon):
sigma14 = atom_type.sigma_14 * length_conv # in md_unit_system
epsilon14 = atom_type.epsilon_14 * ene_conv # in md_unit_system
# Ensure we don't have sigma = 0
if sigma14 == 0.0:
if (epsilon14 == 0.0):
sigma14 = 1.0 # reset sigma = 1
else:
raise ValueError("For atom type '%s', sigma_14 = 0 but "
"epsilon_14 != 0." % name)
else:
sigma14 = None
epsilon14 = None
attributes = { 'class' : name, 'sigma' : str(sigma), 'epsilon' : str(abs(epsilon)) }
if epsilon14 is not None:
attributes['epsilon14'] = str(abs(epsilon14))
if sigma14 is not None:
attributes['sigma14'] = str(sigma14)
etree.SubElement(xml_force, 'Atom', **attributes)
# write NBFIX records
for (atom_types, value) in iteritems(self.nbfix_types):
emin = value[0] * ene_conv
rmin = value[1] * length_conv
# convert to sigma; note that NBFIX types are not rmin/2 but rmin
sigma = rmin/(2**(1.0/6))
etree.SubElement(xml_force, 'NBFixPair', class1=atom_types[0], class2=atom_types[1], sigma=str(sigma), epsilon=str(emin))
@needs_lxml
def _write_omm_DrudeForce(self, xml_root, skip_types):
# Find all atoms with Drude particles.
drude_atoms = []
for residue in list(self.residues.values())+list(self.patches.values()):
for atom in residue.atoms:
if isinstance(atom, DrudeAtom):
drude_atoms.append((atom, residue))
if len(drude_atoms) == 0:
return
if not self.unique_atom_types:
raise ValueError('Drude particles require unique_atom_types')
xml_force = etree.SubElement(xml_root, 'DrudeForce')
alpha_scale = (1*u.angstrom/u.nanometers)**3
for atom, residue in drude_atoms:
attributes = { 'type1' : self._get_mm_atom_type(atom, residue, True), 'type2' : self._get_mm_atom_type(atom, residue),
'charge' : str(atom.drude_charge), 'polarizability' : str(abs(alpha_scale*atom.alpha)),
'thole' : str(atom.thole) }
if atom.anisotropy is not None:
aniso = atom.anisotropy
attributes['type3'] = self._get_mm_atom_type(aniso.atom2, residue)
attributes['type4'] = self._get_mm_atom_type(aniso.atom3, residue)
attributes['type5'] = self._get_mm_atom_type(aniso.atom4, residue)
attributes['aniso12'] = str(aniso.a11)
attributes['aniso34'] = str(aniso.a22)
etree.SubElement(xml_force, 'Particle', **attributes)
def _write_omm_scripts(self, dest, skip_types):
# Not currently implemented, so throw an exception if any unsupported
# options are specified
if self.combining_rule == 'geometric':
raise NotImplementedError('Geometric combining rule not currently supported.')