"""Foyer ForceField class and utility methods."""
import collections
import glob
import itertools
import os
import re
import warnings
import xml.etree.ElementTree as ET
from copy import copy
from tempfile import NamedTemporaryFile
from typing import Callable, Iterable, List
import numpy as np
import openmm as mm
import openmm.app.element as elem
import openmm.unit as u
import parmed as pmd
from openmm import app
from openmm.app.forcefield import (
AllBonds,
CutoffNonPeriodic,
HAngles,
HarmonicAngleGenerator,
HarmonicBondGenerator,
HBonds,
NoCutoff,
NonbondedGenerator,
PeriodicTorsionGenerator,
RBTorsionGenerator,
_convertParameterToNumber,
)
from parmed.gromacs.gromacstop import _Defaults
from pkg_resources import iter_entry_points, resource_filename
import foyer.element as custom_elem
from foyer import smarts
from foyer.atomtyper import find_atomtypes
from foyer.exceptions import (
FoyerError,
MissingForceError,
MissingParametersError,
UnimplementedCombinationRuleError,
)
from foyer.utils.external import get_ref
from foyer.utils.io import has_mbuild, import_
from foyer.utils.misc import validate_type
from foyer.validator import Validator
from foyer.xml_writer import write_foyer
force_for = {
NonbondedGenerator: "NonBondedForce",
HarmonicBondGenerator: "HarmonicBondForce",
HarmonicAngleGenerator: "HarmonicAngleForce",
PeriodicTorsionGenerator: "PeriodicTorsionForce",
RBTorsionGenerator: "RBTorsionForce",
}
def preprocess_forcefield_files(forcefield_files=None):
"""Pre-process foyer Forcefield XML files."""
if forcefield_files is None:
return None
preprocessed_files = []
for xml_file in forcefield_files:
if not hasattr(xml_file, "read"):
f = open(xml_file)
_, suffix = os.path.split(xml_file)
else:
f = xml_file
suffix = ""
# read and preprocess
xml_contents = f.read()
f.close()
xml_contents = re.sub(
r"(def\w*=\w*[\"\'])(.*)([\"\'])",
lambda m: m.group(1)
+ re.sub(r"&(?!amp;)", r"&", m.group(2))
+ m.group(3),
xml_contents,
)
try:
"""
Sort topology objects by precedence, defined by the number of
`type` attributes specified, where a `type` attribute indicates
increased specificity as opposed to use of `class`
"""
root = ET.fromstring(xml_contents)
for element in root:
if "Force" in element.tag:
element[:] = sorted(
element,
key=lambda child: (
-1
* len(
[
attr_name
for attr_name in child.keys()
if "type" in attr_name
]
)
),
)
xml_contents = ET.tostring(root, method="xml").decode()
except ET.ParseError:
"""
Provide the user with a warning if sorting could not be performed.
This indicates a bad XML file, which will be passed on to the
Validator to yield a more descriptive error message.
"""
warnings.warn(
"Invalid XML detected. Could not auto-sort topology "
"objects by precedence."
)
# write to temp file
temp_file = NamedTemporaryFile(suffix=suffix, delete=False)
with open(temp_file.name, "w") as temp_f:
temp_f.write(xml_contents)
# append temp file name to list
preprocessed_files.append(temp_file.name)
return preprocessed_files
def get_available_forcefield_loaders() -> List[Callable]:
"""Get a list of available force field loader functions."""
available_ff_paths = []
for entry_point in iter_entry_points(group="foyer.forcefields"):
available_ff_paths.append(entry_point.load())
return available_ff_paths
def generate_topology(non_omm_topology, non_element_types=None, residues=None):
"""Create an OpenMM Topology from another supported topology structure."""
if non_element_types is None:
non_element_types = set()
if isinstance(non_omm_topology, pmd.Structure):
return _topology_from_parmed(non_omm_topology, non_element_types)
elif has_mbuild:
mb = import_("mbuild")
if (non_omm_topology, mb.Compound):
pmd_comp_struct = non_omm_topology.to_parmed(residues=residues)
return _topology_from_parmed(pmd_comp_struct, non_element_types)
else:
raise FoyerError(
"Unknown topology format: {}\n"
"Supported formats are: "
'"parmed.Structure", '
'"mbuild.Compound", '
'"openmm.app.Topology"'.format(non_omm_topology)
)
def _structure_from_residue(residue, parent_structure):
"""Convert a ParmEd Residue to an equivalent Structure."""
structure = pmd.Structure()
orig_to_copy = (
dict()
) # Clone a lot of atoms to avoid any of parmed's tracking
for atom in residue.atoms:
new_atom = copy(atom)
new_atom._idx = atom.idx
orig_to_copy[atom] = new_atom
structure.add_atom(
new_atom, resname=residue.name, resnum=residue.number
)
for bond in parent_structure.bonds:
if bond.atom1 in residue.atoms and bond.atom2 in residue.atoms:
structure.bonds.append(
pmd.Bond(orig_to_copy[bond.atom1], orig_to_copy[bond.atom2])
)
return structure
def _topology_from_parmed(structure, non_element_types):
"""Convert a ParmEd Structure to an OpenMM Topology."""
topology = app.Topology()
residues = dict()
for pmd_residue in structure.residues:
chain = topology.addChain()
omm_residue = topology.addResidue(pmd_residue.name, chain)
# Index ParmEd residues on name & number, no other info i.e. chain
residues[(pmd_residue.name, pmd_residue.idx)] = omm_residue
atoms = dict() # pmd.Atom: omm.Atom
for pmd_atom in structure.atoms:
name = pmd_atom.name
if pmd_atom.name in non_element_types:
element = non_element_types[pmd_atom.name]
else:
if (
isinstance(pmd_atom.atomic_number, int)
and pmd_atom.atomic_number != 0
):
element = elem.Element.getByAtomicNumber(pmd_atom.atomic_number)
else:
element = elem.Element.getBySymbol(pmd_atom.name)
omm_atom = topology.addAtom(
name,
element,
residues[(pmd_atom.residue.name, pmd_atom.residue.idx)],
)
omm_atom.id = pmd_atom.id
atoms[pmd_atom] = omm_atom
omm_atom.bond_partners = []
for bond in structure.bonds:
atom1 = atoms[bond.atom1]
atom2 = atoms[bond.atom2]
topology.addBond(atom1, atom2)
atom1.bond_partners.append(atom2)
atom2.bond_partners.append(atom1)
if structure.box_vectors and np.any(
[x._value for x in structure.box_vectors]
):
topology.setPeriodicBoxVectors(structure.box_vectors)
positions = structure.positions
return topology, positions
def _topology_from_residue(res):
"""Convert an openmm.app.Topology.Residue to openmm.app.Topology.
Parameters
----------
res : openmm.app.Topology.Residue
An individual residue in an openmm.app.Topology
Returns
-------
topology : openmm.app.Topology
The generated topology
"""
topology = app.Topology()
chain = topology.addChain()
new_res = topology.addResidue(res.name, chain)
atoms = dict() # { omm.Atom in res : omm.Atom in *new* topology }
for res_atom in res.atoms():
topology_atom = topology.addAtom(
name=res_atom.name, element=res_atom.element, residue=new_res
)
atoms[res_atom] = topology_atom
topology_atom.bond_partners = []
for bond in res.bonds():
atom1 = atoms[bond.atom1]
atom2 = atoms[bond.atom2]
topology.addBond(atom1, atom2)
atom1.bond_partners.append(atom2)
atom2.bond_partners.append(atom1)
return topology
def _check_independent_residues(structure):
"""Check to see if residues will constitute independent graphs."""
for res in structure.residues:
atoms_in_residue = set([*res.atoms])
bond_partners_in_residue = [
item
for sublist in [atom.bond_partners for atom in res.atoms]
for item in sublist
]
# Handle the case of a 'residue' with no neighbors
if not bond_partners_in_residue:
continue
if set(atoms_in_residue) != set(bond_partners_in_residue):
return False
return True
def _unwrap_typemap(structure, residue_map):
master_typemap = {
atom.idx: {"whitelist": set(), "blacklist": set(), "atomtype": None}
for atom in structure.atoms
}
for res in structure.residues:
for res_ref, val in residue_map.items():
if id(res.name) == id(res_ref):
for i, atom in enumerate(res.atoms):
master_typemap[int(atom.idx)]["atomtype"] = val[i][
"atomtype"
]
return master_typemap
def _separate_urey_bradleys(system, topology):
"""Separate urey bradley bonds from harmonic bonds in OpenMM System.
Parameters
----------
topology : openmm.app.Topology
Molecular structure to find atom types of
system : openmm System
"""
atoms = [a for a in topology.atoms()]
bonds = [b for b in topology.bonds()]
ub_force = mm.HarmonicBondForce()
harmonic_bond_force = mm.HarmonicBondForce()
for force_idx, force in enumerate(system.getForces()):
if isinstance(force, mm.HarmonicBondForce):
for bond_idx in range(force.getNumBonds()):
if (
atoms[force.getBondParameters(bond_idx)[0]],
atoms[force.getBondParameters(bond_idx)[1]],
) not in bonds and (
atoms[force.getBondParameters(bond_idx)[1]],
atoms[force.getBondParameters(bond_idx)[0]],
) not in bonds:
ub_force.addBond(*force.getBondParameters(bond_idx))
else:
harmonic_bond_force.addBond(
*force.getBondParameters(bond_idx)
)
system.removeForce(force_idx)
system.addForce(harmonic_bond_force)
system.addForce(ub_force)
def _error_or_warn(error, msg):
"""Raise an error or warning if topology objects are not fully parameterized.
Parameters
----------
error : bool
If True, raise an error, else raise a warning
msg : str
The message to be provided with the error or warning
"""
if error:
raise Exception(msg)
else:
warnings.warn(msg)
def _check_bonds(data, structure, verbose, assert_bond_params):
"""Check if any bonds lack paramters."""
if data.bonds:
missing = [b for b in structure.bonds if b.type is None]
if missing:
if verbose:
for bond in missing:
a1 = bond.atom1
a2 = bond.atom2
print(
f"Missing bond with ids {(a1.idx, a2.idx)} and types "
f"{[a1.type, a2.type]}."
)
nmissing = len(structure.bonds) - len(missing)
msg = (
"Parameters have not been assigned to all bonds. "
"Total system bonds: {}, Parametrized bonds: {}"
"".format(len(structure.bonds), nmissing)
)
_error_or_warn(assert_bond_params, msg)
def _check_angles(data, structure, verbose, assert_angle_params):
"""Check if all angles were found and parametrized."""
if verbose:
for omm_ids in data.angles:
missing_angle = True
for pmd_angle in structure.angles:
pmd_ids = (
pmd_angle.atom1.idx,
pmd_angle.atom2.idx,
pmd_angle.atom3.idx,
)
if pmd_ids == omm_ids:
missing_angle = False
if missing_angle:
print(
"Missing angle with ids {} and types {}.".format(
omm_ids, [structure.atoms[idx].type for idx in omm_ids]
)
)
if data.angles and (len(data.angles) != len(structure.angles)):
msg = (
"Parameters have not been assigned to all angles. Total "
"system angles: {}, Parameterized angles: {}"
"".format(len(data.angles), len(structure.angles))
)
_error_or_warn(assert_angle_params, msg)
def _check_dihedrals(
data, structure, verbose, assert_dihedral_params, assert_improper_params
):
"""Check if all dihedrals, including impropers, were found and parametrized."""
proper_dihedrals = [
dihedral for dihedral in structure.dihedrals if not dihedral.improper
]
if verbose:
for omm_ids in data.propers:
missing_dihedral = True
for pmd_proper in proper_dihedrals:
pmd_ids = (
pmd_proper.atom1.idx,
pmd_proper.atom2.idx,
pmd_proper.atom3.idx,
pmd_proper.atom4.idx,
)
if pmd_ids == omm_ids:
missing_dihedral = False
for pmd_proper in structure.rb_torsions:
pmd_ids = (
pmd_proper.atom1.idx,
pmd_proper.atom2.idx,
pmd_proper.atom3.idx,
pmd_proper.atom4.idx,
)
if pmd_ids == omm_ids:
missing_dihedral = False
if missing_dihedral:
print(
"Missing dihedral with ids {} and types {}.".format(
omm_ids, [structure.atoms[idx].type for idx in omm_ids]
)
)
if data.propers and len(data.propers) != len(proper_dihedrals) + len(
structure.rb_torsions
):
if data.propers and len(data.propers) < len(proper_dihedrals) + len(
structure.rb_torsions
):
msg = (
"Parameters have been assigned to all proper dihedrals. "
"However, there are more parameterized dihedrals ({}) "
"than total system dihedrals ({}). "
"This may be due to having multiple periodic dihedrals "
"for a single system dihedral.".format(
len(proper_dihedrals) + len(structure.rb_torsions),
len(data.propers),
)
)
warnings.warn(msg)
else:
msg = (
"Parameters have not been assigned to all proper dihedrals. "
"Total system dihedrals: {}, Parameterized dihedrals: {}. "
"Note that if your system contains torsions of Ryckaert-"
"Bellemans functional form, all of these torsions are "
"processed as propers.".format(
len(data.propers),
len(proper_dihedrals) + len(structure.rb_torsions),
)
)
_error_or_warn(assert_dihedral_params, msg)
improper_dihedrals = [
dihedral for dihedral in structure.dihedrals if dihedral.improper
]
if data.impropers and len(data.impropers) != len(improper_dihedrals) + len(
structure.impropers
):
msg = (
"Parameters have not been assigned to all impropers. Total "
"system impropers: {}, Parameterized impropers: {}. "
"Note that if your system contains torsions of Ryckaert-"
"Bellemans functional form, all of these torsions are "
"processed as propers".format(
len(data.impropers),
len(improper_dihedrals) + len(structure.impropers),
)
)
_error_or_warn(assert_improper_params, msg)
[docs]class Forcefield(app.ForceField):
"""Specialization of OpenMM's Forcefield allowing SMARTS based atomtyping.
Parameters
----------
forcefield_files : list of str, optional, default=None
List of forcefield files to load.
name : str, optional, default=None
Name of a forcefield to load that is packaged within foyer.
"""
def __init__(
self, forcefield_files=None, name=None, validation=True, debug=False
):
self.atomTypeDefinitions = dict()
self.atomTypeOverrides = dict()
self.atomTypeDesc = dict()
self.atomTypeRefs = dict()
self.atomTypeClasses = dict()
self.atomTypeElements = dict()
self._included_forcefields = dict()
self.non_element_types = dict()
self._version = None
self._name = None
self._combining_rule = None
all_files_to_load = []
if forcefield_files is not None:
if isinstance(forcefield_files, (list, tuple, set)):
for file in forcefield_files:
all_files_to_load.append(file)
else:
all_files_to_load.append(forcefield_files)
if name is not None:
try:
file = self.included_forcefields[name]
except KeyError:
raise IOError("Forcefield {} cannot be found".format(name))
else:
all_files_to_load.append(file)
preprocessed_files = preprocess_forcefield_files(all_files_to_load)
if validation:
for ff_file_name in preprocessed_files:
Validator(ff_file_name, debug)
super(Forcefield, self).__init__(*preprocessed_files)
if len(preprocessed_files) == 1:
self._version = self._parse_version_number(preprocessed_files[0])
self._name = self._parse_name(preprocessed_files[0])
self._combining_rule = self._parse_combining_rule(
preprocessed_files[0]
)
elif len(preprocessed_files) > 1:
self._version = [
self._parse_version_number(f) for f in preprocessed_files
]
self._name = [self._parse_name(f) for f in preprocessed_files]
self._combining_rule = [
self._parse_combining_rule(f) for f in preprocessed_files
]
if len(set(self._combining_rule)) == 1:
self._combining_rule = self._combining_rule[0]
else:
raise FoyerError(
"Inconsistent combining_rule among loaded forcefield files"
)
for fp in preprocessed_files:
os.remove(fp)
self.parser = smarts.SMARTS(self.non_element_types)
self._system_data = None
@property
def version(self):
"""Return version number of the force field XML file."""
return self._version
@property
def name(self):
"""Return the name of the force field XML."""
return self._name
@property
def combining_rule(self):
"""Return the combining rule of this force field."""
return self._combining_rule
@property
def included_forcefields(self):
"""Return a dictionary mappying for all included forcefields."""
if any(self._included_forcefields):
return self._included_forcefields
ff_dir = resource_filename("foyer", "forcefields")
ff_filepaths = set(glob.glob(os.path.join(ff_dir, "xml/*.xml")))
for ff_filepath in ff_filepaths:
_, ff_file = os.path.split(ff_filepath)
basename, _ = os.path.splitext(ff_file)
self._included_forcefields[basename] = ff_filepath
return self._included_forcefields
@property
def lj14scale(self):
"""Get LJ 1-4 scale for this forcefield."""
try:
non_bonded_force_gen = self.get_generator(
ff=self, gen_type=NonbondedGenerator
)
except MissingForceError:
raise AttributeError(
"Cannot get lj14Scale for the forcefield "
"because it doesn't have NonBondedForce."
)
return non_bonded_force_gen.lj14scale
@property
def coulomb14scale(self):
"""Get Coulombic 1-4 scale for this forcefield."""
try:
non_bonded_force_gen = self.get_generator(
ff=self, gen_type=NonbondedGenerator
)
except MissingForceError:
raise AttributeError(
"Cannot get coulomb14scale for the Forcefield "
"because it doesn't have NonBondedForce."
)
return non_bonded_force_gen.coulomb14scale
def _parse_version_number(self, forcefield_file):
with open(forcefield_file, "r") as f:
tree = ET.parse(f)
root = tree.getroot()
try:
return root.attrib["version"]
except KeyError:
warnings.warn(
"No force field version number found in force field XML file."
)
return None
def _parse_name(self, forcefield_file):
with open(forcefield_file, "r") as f:
tree = ET.parse(f)
root = tree.getroot()
try:
return root.attrib["name"]
except KeyError:
warnings.warn(
"No force field name found in force field XML file."
)
return None
def _parse_combining_rule(self, forcefield_file):
with open(forcefield_file, "r") as f:
tree = ET.parse(f)
root = tree.getroot()
try:
return root.attrib["combining_rule"]
except KeyError:
warnings.warn(
"No combining rule found in force field XML file."
)
return "lorentz"
def _create_element(self, element, mass):
if not isinstance(element, elem.Element):
try:
element = elem.get_by_symbol(element)
except KeyError:
# Enables support for non-atomistic "element types"
if element not in self.non_element_types:
warnings.warn(
"Non-atomistic element type detected. "
"Creating custom element for {}".format(element)
)
element = custom_elem.Element(
number=0, mass=mass, name=element, symbol=element
)
else:
return element, False
return element, True
[docs] def registerAtomType(self, parameters):
"""Register a new atom type."""
name = parameters["name"]
if name in self._atomTypes:
raise ValueError(
"Found multiple definitions for atom type: " + name
)
atom_class = parameters["class"]
mass = _convertParameterToNumber(parameters["mass"])
element = None
if "element" in parameters:
element, custom = self._create_element(parameters["element"], mass)
if custom:
self.non_element_types[element.symbol] = element
self._atomTypes[name] = self.__class__._AtomType(
name, atom_class, mass, element
)
if atom_class in self._atomClasses:
type_set = self._atomClasses[atom_class]
else:
type_set = set()
self._atomClasses[atom_class] = type_set
type_set.add(name)
self._atomClasses[""].add(name)
name = parameters["name"]
if "def" in parameters:
self.atomTypeDefinitions[name] = parameters["def"]
if "overrides" in parameters:
overrides = set(
atype.strip() for atype in parameters["overrides"].split(",")
)
if overrides:
self.atomTypeOverrides[name] = overrides
if "desc" in parameters:
self.atomTypeDesc[name] = parameters["desc"]
if "doi" in parameters:
dois = set(doi.strip() for doi in parameters["doi"].split(","))
self.atomTypeRefs[name] = dois
if "element" in parameters:
self.atomTypeElements[name] = parameters["element"]
if "class" in parameters:
self.atomTypeClasses[name] = parameters["class"]
[docs] def apply(
self,
structure,
references_file=None,
use_residue_map=True,
assert_bond_params=True,
assert_angle_params=True,
assert_dihedral_params=True,
assert_improper_params=False,
verbose=False,
*args,
**kwargs,
):
"""Apply the force field to a molecular structure.
Parameters
----------
structure : parmed.Structure or mbuild.Compound
Molecular structure to apply the force field to
references_file : str, optional, default=None
Name of file where force field references will be written (in Bibtex
format)
use_residue_map : boolean, optional, default=True
Store atomtyped topologies of residues to a dictionary that maps
them to residue names. Each topology, including atomtypes, will be
copied to other residues with the same name. This avoids repeatedly
calling the subgraph isomorphism on idential residues and should
result in better performance for systems with many identical
residues, i.e. a box of water. Note that for this to be applied to
independent molecules, they must each be saved as different
residues in the topology.
assert_bond_params : bool, optional, default=True
If True, Foyer will exit if parameters are not found for all system
bonds.
assert_angle_params : bool, optional, default=True
If True, Foyer will exit if parameters are not found for all system
angles.
assert_dihedral_params : bool, optional, default=True
If True, Foyer will exit if parameters are not found for all system
proper dihedrals.
assert_improper_params : bool, optional, default=False
If True, Foyer will exit if parameters are not found for all system
improper dihedrals.
verbose : bool, optional, default=False
If True, Foyer will print debug-level information about notable or
potentially problematic details it encounters.
"""
if self.atomTypeDefinitions == {}:
raise FoyerError(
"Attempting to atom-type using a force field "
"with no atom type defitions."
)
if not isinstance(structure, pmd.Structure):
mb = import_("mbuild")
if isinstance(structure, mb.Compound):
structure = structure.to_parmed(**kwargs)
typemap = self.run_atomtyping(
structure, use_residue_map=use_residue_map, **kwargs
)
self._apply_typemap(structure, typemap)
structure = self.parametrize_system(
structure=structure,
references_file=references_file,
assert_bond_params=assert_bond_params,
assert_angle_params=assert_angle_params,
assert_dihedral_params=assert_dihedral_params,
assert_improper_params=assert_improper_params,
verbose=verbose,
*args,
**kwargs,
)
# Start storing scaling factors and combining rule to parmed structure
# Utilizing parmed gromactop's _Default class
# Note: nb_func = 1 (LJ) or 2 (Buckingham), foyer forcefield only support 1 at the moment
# Combining rule follow GROMACS scheme, where "lorentz" is 2 and "geometric" is 3
combining_rules = {"lorentz": 2, "geometric": 3}
gen_pairs = "yes" if structure.adjusts else "no"
try:
lj14scale = self.lj14scale
coulomb14scale = self.coulomb14scale
structure.defaults = _Defaults(
nbfunc=1,
comb_rule=combining_rules[self.combining_rule],
gen_pairs=gen_pairs,
fudgeLJ=lj14scale,
fudgeQQ=coulomb14scale,
)
except AttributeError:
warnings.warn(
"Missing lj14scale or coulomb14scale, could not set structure metadata."
)
structure.defaults = None
return structure
[docs] def run_atomtyping(self, structure, use_residue_map=True, **kwargs):
"""Atomtype the topology.
Parameters
----------
structure : parmed.structure.Structure
Molecular structure to find atom types of
use_residue_map : boolean, optional, default=True
Store atomtyped topologies of residues to a dictionary that maps
them to residue names. Each topology, including atomtypes, will be
copied to other residues with the same name. This avoids repeatedly
calling the subgraph isomorphism on idential residues and should
result in better performance for systems with many identical
residues, i.e. a box of water. Note that for this to be applied to
independent molecules, they must each be saved as different
residues in the topology.
"""
if use_residue_map:
independent_residues = _check_independent_residues(structure)
if independent_residues:
residue_map = dict()
# Need to call this only once and store results for later id() comparisons
for res_id, res in enumerate(structure.residues):
if (
structure.residues[res_id].name
not in residue_map.keys()
):
tmp_res = _structure_from_residue(res, structure)
typemap = find_atomtypes(tmp_res, forcefield=self)
residue_map[res.name] = typemap
typemap = _unwrap_typemap(structure, residue_map)
else:
typemap = find_atomtypes(structure, forcefield=self)
else:
typemap = find_atomtypes(structure, forcefield=self)
return typemap
[docs] def parametrize_system(
self,
structure=None,
references_file=None,
assert_bond_params=True,
assert_angle_params=True,
assert_dihedral_params=True,
assert_improper_params=False,
verbose=False,
*args,
**kwargs,
):
"""Create system based on resulting typemapping."""
topology, positions = _topology_from_parmed(
structure, self.non_element_types
)
system = self.createSystem(topology, *args, **kwargs)
_separate_urey_bradleys(system, topology)
data = self._system_data
structure = pmd.openmm.load_topology(topology=topology, system=system)
structure.bonds.sort(key=lambda x: x.atom1.idx)
structure.positions = positions
box_vectors = topology.getPeriodicBoxVectors()
if box_vectors is not None:
structure.box_vectors = box_vectors
_check_bonds(data, structure, verbose, assert_bond_params)
_check_angles(data, structure, verbose, assert_angle_params)
_check_dihedrals(
data,
structure,
verbose,
assert_dihedral_params,
assert_improper_params,
)
if references_file:
atom_types = set(atom.type for atom in structure.atoms)
self._write_references_to_file(atom_types, references_file)
try:
structure.combining_rule = self.combining_rule
except ValueError as e:
# Parmed raises ValueError: combining_rule must be 'lorentz' or 'geometric'
if str(e).startswith("combining_rule must be"):
raise UnimplementedCombinationRuleError(
f"Combination rule {self.combining_rule} is not implemented"
)
if self.combining_rule == "geometric":
self._patch_parmed_adjusts(
structure, combining_rule=self.combining_rule
)
total_charge = sum([atom.charge for atom in structure.atoms])
if not np.allclose(total_charge, 0):
warnings.warn(
"Parametrized structure has non-zero charge."
"Structure's total charge: {}".format(total_charge)
)
return structure
[docs] def createSystem(
self,
topology,
nonbondedMethod=NoCutoff,
nonbondedCutoff=1.0 * u.nanometer,
constraints=None,
rigidWater=True,
removeCMMotion=True,
hydrogenMass=None,
switchDistance=None,
**args,
):
"""Construct an OpenMM System representing a Topology with this force field.
Parameters
----------
topology : Topology
The Topology for which to create a System
nonbondedMethod : object=NoCutoff
The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
nonbondedCutoff : distance=1*nanometer
The cutoff distance to use for nonbonded interactions
constraints : object=None
Specifies which bonds and angles should be implemented with constraints.
Allowed values are None, HBonds, AllBonds, or HAngles.
rigidWater : boolean=True
If true, water molecules will be fully rigid regardless of the value
passed for the constraints argument
removeCMMotion : boolean=True
If true, a CMMotionRemover will be added to the System
hydrogenMass : mass=None
The mass to use for hydrogen atoms bound to heavy atoms. Any mass
added to a hydrogen is subtracted from the heavy atom to keep
their total mass the same.
switchDistance : float=None
The distance at which the potential energy switching function is turned on for
args
Arbitrary additional keyword arguments may also be specified.
This allows extra parameters to be specified that are specific to
particular force fields.
Returns
-------
system
the newly created System
"""
args["switchDistance"] = switchDistance
# Overwrite previous _SystemData object
data = app.ForceField._SystemData(topology)
self._system_data = data
# TODO: Better way to lookup nonbonded parameters...?
nonbonded_params = None
for generator in self.getGenerators():
if isinstance(generator, NonbondedGenerator):
nonbonded_params = generator.params.paramsForType
break
for chain in topology.chains():
for res in chain.residues():
for atom in res.atoms():
data.atomType[atom] = atom.id
if nonbonded_params:
params = nonbonded_params[atom.id]
data.atomParameters[atom] = params
# Create the System and add atoms
sys = mm.System()
for atom in topology.atoms():
# Look up the atom type name, returning a helpful error message if it cannot be found.
if atom not in data.atomType:
raise Exception(
"Could not identify atom type for atom '%s'." % str(atom)
)
typename = data.atomType[atom]
# Look up the type name in the list of registered atom types, returning a helpful error message if it cannot be found.
if typename not in self._atomTypes:
msg = (
"Could not find typename '%s' for atom '%s' in list of known atom types.\n"
% (typename, str(atom))
)
msg += "Known atom types are: %s" % str(self._atomTypes.keys())
raise Exception(msg)
# Add the particle to the OpenMM system.
mass = self._atomTypes[typename].mass
sys.addParticle(mass)
# Adjust hydrogen masses if requested.
if hydrogenMass is not None:
if not u.is_quantity(hydrogenMass):
hydrogenMass *= u.dalton
for atom1, atom2 in topology.bonds():
if atom1.element == elem.hydrogen:
(atom1, atom2) = (atom2, atom1)
if atom2.element == elem.hydrogen and atom1.element not in (
elem.hydrogen,
None,
):
transfer_mass = hydrogenMass - sys.getParticleMass(
atom2.index
)
sys.setParticleMass(atom2.index, hydrogenMass)
mass = sys.getParticleMass(atom1.index) - transfer_mass
sys.setParticleMass(atom1.index, mass)
# Set periodic boundary conditions.
box_vectors = topology.getPeriodicBoxVectors()
if box_vectors is not None:
sys.setDefaultPeriodicBoxVectors(
box_vectors[0], box_vectors[1], box_vectors[2]
)
elif nonbondedMethod not in [NoCutoff, CutoffNonPeriodic]:
raise ValueError(
"Requested periodic boundary conditions for a "
"Topology that does not specify periodic box "
"dimensions"
)
# Make a list of all unique angles
unique_angles = set()
for bond in data.bonds:
for atom in data.bondedToAtom[bond.atom1]:
if atom != bond.atom2:
if atom < bond.atom2:
unique_angles.add((atom, bond.atom1, bond.atom2))
else:
unique_angles.add((bond.atom2, bond.atom1, atom))
for atom in data.bondedToAtom[bond.atom2]:
if atom != bond.atom1:
if atom > bond.atom1:
unique_angles.add((bond.atom1, bond.atom2, atom))
else:
unique_angles.add((atom, bond.atom2, bond.atom1))
data.angles = sorted(list(unique_angles))
# Make a list of all unique proper torsions
unique_propers = set()
for angle in data.angles:
for atom in data.bondedToAtom[angle[0]]:
if atom not in angle:
if atom < angle[2]:
unique_propers.add((atom, angle[0], angle[1], angle[2]))
else:
unique_propers.add((angle[2], angle[1], angle[0], atom))
for atom in data.bondedToAtom[angle[2]]:
if atom not in angle:
if atom > angle[0]:
unique_propers.add((angle[0], angle[1], angle[2], atom))
else:
unique_propers.add((atom, angle[2], angle[1], angle[0]))
data.propers = sorted(list(unique_propers))
# Make a list of all unique improper torsions
for atom in range(len(data.bondedToAtom)):
bonded_to = data.bondedToAtom[atom]
if len(bonded_to) > 2:
for subset in itertools.combinations(bonded_to, 3):
data.impropers.append(
(atom, subset[0], subset[1], subset[2])
)
# Identify bonds that should be implemented with constraints
if constraints == AllBonds or constraints == HAngles:
for bond in data.bonds:
bond.isConstrained = True
elif constraints == HBonds:
for bond in data.bonds:
atom1 = data.atoms[bond.atom1]
atom2 = data.atoms[bond.atom2]
bond.isConstrained = atom1.name.startswith(
"H"
) or atom2.name.startswith("H")
if rigidWater:
for bond in data.bonds:
atom1 = data.atoms[bond.atom1]
atom2 = data.atoms[bond.atom2]
if atom1.residue.name == "HOH" and atom2.residue.name == "HOH":
bond.isConstrained = True
# Identify angles that should be implemented with constraints
if constraints == HAngles:
for angle in data.angles:
atom1 = data.atoms[angle[0]]
atom2 = data.atoms[angle[1]]
atom3 = data.atoms[angle[2]]
numH = 0
if atom1.name.startswith("H"):
numH += 1
if atom3.name.startswith("H"):
numH += 1
data.isAngleConstrained.append(
numH == 2 or (numH == 1 and atom2.name.startswith("O"))
)
else:
data.isAngleConstrained = len(data.angles) * [False]
if rigidWater:
for i in range(len(data.angles)):
angle = data.angles[i]
atom1 = data.atoms[angle[0]]
atom2 = data.atoms[angle[1]]
atom3 = data.atoms[angle[2]]
if (
atom1.residue.name == "HOH"
and atom2.residue.name == "HOH"
and atom3.residue.name == "HOH"
):
data.isAngleConstrained[i] = True
# Add virtual sites
for atom in data.virtualSites:
(site, atoms, excludeWith) = data.virtualSites[atom]
index = atom.index
data.excludeAtomWith[excludeWith].append(index)
if site.type == "average2":
sys.setVirtualSite(
index,
mm.TwoParticleAverageSite(
atoms[0], atoms[1], site.weights[0], site.weights[1]
),
)
elif site.type == "average3":
sys.setVirtualSite(
index,
mm.ThreeParticleAverageSite(
atoms[0],
atoms[1],
atoms[2],
site.weights[0],
site.weights[1],
site.weights[2],
),
)
elif site.type == "outOfPlane":
sys.setVirtualSite(
index,
mm.OutOfPlaneSite(
atoms[0],
atoms[1],
atoms[2],
site.weights[0],
site.weights[1],
site.weights[2],
),
)
elif site.type == "localCoords":
local_coord_site = mm.LocalCoordinatesSite(
atoms[0],
atoms[1],
atoms[2],
mm.Vec3(
site.originWeights[0],
site.originWeights[1],
site.originWeights[2],
),
mm.Vec3(
site.xWeights[0], site.xWeights[1], site.xWeights[2]
),
mm.Vec3(
site.yWeights[0], site.yWeights[1], site.yWeights[2]
),
mm.Vec3(
site.localPos[0], site.localPos[1], site.localPos[2]
),
)
sys.setVirtualSite(index, local_coord_site)
# Add forces to the System
for force in self._forces:
force.createForce(sys, data, nonbondedMethod, nonbondedCutoff, args)
if removeCMMotion:
sys.addForce(mm.CMMotionRemover())
# Let force generators do postprocessing
for force in self._forces:
if "postprocessSystem" in dir(force):
force.postprocessSystem(sys, data, args)
# Execute scripts found in the XML files.
for script in self._scripts:
exec(script, locals())
return sys
def _apply_typemap(self, structure, typemap):
"""Add atomtypes to the topology according to the typemap."""
for atom in structure.atoms:
atom.id = typemap[atom.idx]["atomtype"]
if not all([a.id for a in structure.atoms]):
raise ValueError("Not all atoms in topology have atom types")
def _prepare_structure(self, topology, **kwargs):
"""Separate positions and other topological information."""
if not isinstance(topology, app.Topology):
residues = kwargs.get("residues")
topology, positions = generate_topology(
topology, self.non_element_types, residues=residues
)
else:
positions = np.empty(shape=(topology.getNumAtoms(), 3))
positions[:] = np.nan
return topology, positions
def _patch_parmed_adjusts(self, structure, combining_rule="geometric"):
if combining_rule != "geometric":
raise UnimplementedCombinationRuleError()
for adj in structure.adjusts:
sig1 = adj.atom1.sigma
sig2 = adj.atom2.sigma
sig_geo = (sig1 * sig2) ** 0.5
adj.type.sigma = sig_geo
return structure
def _write_references_to_file(self, atom_types, references_file):
atomtype_references = {}
for atype in atom_types:
try:
atomtype_references[atype] = self.atomTypeRefs[atype]
except KeyError:
warnings.warn(
"Reference not found for atom type '{}'." "".format(atype)
)
unique_references = collections.defaultdict(list)
for atomtype, dois in atomtype_references.items():
for doi in dois:
unique_references[doi].append(atomtype)
unique_references = collections.OrderedDict(
sorted(unique_references.items())
)
with open(references_file, "w") as f:
for doi, atomtypes in unique_references.items():
url = "http://api.crossref.org/works/{}/transform/application/x-bibtex".format(
doi
)
headers = {"accept": "application/x-bibtex"}
bibtex_ref = get_ref(url, headers=headers)
if bibtex_ref is None:
warnings.warn("Could not get ref for doi {}".format(doi))
continue
else:
bibtex_text = bibtex_ref.text
note = (
",\n\tnote = {Parameters for atom types: "
+ ", ".join(sorted(atomtypes))
+ "}"
)
bibtex_text = bibtex_text[:-2] + note + bibtex_text[-2:]
f.write("{}\n".format(bibtex_text))
[docs] def get_parameters(self, group, key, keys_are_atom_classes=False):
"""Get parameters for a specific group of Forces in this Forcefield.
Parameters
----------
group: str
One of {"atoms", "harmonic_bonds", "harmonic_angles", "periodic_propers", "periodic_impropers", "rb_propers", "rb_impropers"}.
Note that these entries are case insensitive
key: str, list of str
The atom type(s)/class(es) to extract parameters for
keys_are_atom_classes: bool, default=False
If True, the entries in key are considered to be atom classes rather than atom types
Examples
--------
Following shows some example usage of this function:
>>> from foyer import Forcefield
>>> ff = Forcefield(name='oplsaa')
>>> ff.get_parameters('atoms', key='opls_235') # NonBonded Parameters
{'charge': -0.18, 'sigma': 0.35, 'epsilon': 0.276144}
>>> ff.get_parameters('rb_propers', ['opls_137', 'opls_209', 'opls_193', 'opls_154']) # RBPropers
{'c0': 2.87441, 'c1': 0.58158, 'c2': 2.092, 'c3': -5.54799, 'c4': 0.0, 'c5': 0.0}
Returns
-------
dict or list of dict
A dictionary (or an item in the list of dictionaries) with parameter names as key and the parameter value as value as it pertains to
specific force. Note that the keys of the dictionary can be different for different forces.
Raises
------
MissingParametersError
Raised when parameters are missing from the forcefield or no matching parameters found
MissingForceError
Raised when a particular force generator is missing from the Forcefield
"""
group = group.lower()
param_extractors = {
"atoms": self._extract_non_bonded_params,
"harmonic_bonds": self._extract_harmonic_bond_params,
"harmonic_angles": self._extract_harmonic_angle_params,
"periodic_propers": self._extract_periodic_proper_params,
"periodic_impropers": self._extract_periodic_improper_params,
"rb_propers": self._extract_rb_proper_params,
"rb_impropers": self._extract_rb_improper_params,
}
if group not in param_extractors:
raise ValueError(f"Cannot extract parameters for {group}")
key = (
[key]
if isinstance(key, str) or not isinstance(key, Iterable)
else key
)
validate_type(key, str)
if keys_are_atom_classes:
atom_types = self.map_atom_classes_to_types(key)
try:
return param_extractors[group](atom_types)
except MissingParametersError:
raise MissingParametersError(
f"Missing {group} parameters between between class(es) {key}"
)
return param_extractors[group](key)
def _extract_non_bonded_params(self, atom_type):
"""Return parameters for a specific atom type."""
if isinstance(atom_type, list):
assert len(atom_type) == 1, (
f"NonBondedForce parameters can only be extracted "
f"for one atom type. Provided {len(atom_type)}"
)
atom_type = atom_type[0]
non_bonded_forces_gen = self.get_generator(
ff=self, gen_type=NonbondedGenerator
)
non_bonded_params = non_bonded_forces_gen.params.paramsForType
try:
return non_bonded_params[atom_type]
except KeyError:
raise MissingParametersError(
f"Missing parameters for atom {atom_type}"
)
def _extract_harmonic_bond_params(self, atom_types):
"""Return parameters for a specific HarmonicBondForce between atom types."""
if len(atom_types) != 2:
raise ValueError(
f"HarmonicBond parameters can only "
f"be extracted for two atoms. Provided {len(atom_types)}"
)
harmonic_bond_force_gen = self.get_generator(
ff=self, gen_type=HarmonicBondGenerator
)
atom_1_type = atom_types[0]
atom_2_type = atom_types[1]
bonds_for_atom_type = harmonic_bond_force_gen.bondsForAtomType
params = None
for i in bonds_for_atom_type[atom_1_type]:
types1 = harmonic_bond_force_gen.types1[i]
types2 = harmonic_bond_force_gen.types2[i]
if (atom_1_type in types1 and atom_2_type in types2) or (
atom_2_type in types1 and atom_1_type in types2
):
params = {
"k": harmonic_bond_force_gen.k[i],
"length": harmonic_bond_force_gen.length[i],
}
break
if params is None:
raise MissingParametersError(
f"Missing HarmonicBond parameters for bond between "
f"atoms {atom_1_type} and {atom_2_type}"
)
return params
def _extract_harmonic_angle_params(self, atom_types):
"""Return parameters for a HarmonicAngle between atom types."""
if len(atom_types) != 3:
raise ValueError(
f"HarmonicAngle parameters can only "
f"be extracted for three atoms. Provided {len(atom_types)}"
)
harmonic_angle_force_gen = self.get_generator(
ff=self, gen_type=HarmonicAngleGenerator
)
atom_1_type = atom_types[0]
atom_2_type = atom_types[1]
atom_3_type = atom_types[2]
angles_for_atom_2_type = harmonic_angle_force_gen.anglesForAtom2Type
params = None
for i in angles_for_atom_2_type[atom_2_type]:
types1 = harmonic_angle_force_gen.types1[i]
types2 = harmonic_angle_force_gen.types2[i]
types3 = harmonic_angle_force_gen.types3[i]
if (
atom_1_type in types1
and atom_2_type in types2
and atom_3_type in types3
) or (
atom_3_type in types1
and atom_2_type in types2
and atom_1_type in types3
):
params = {
"theta": harmonic_angle_force_gen.angle[i],
"k": harmonic_angle_force_gen.k[i],
}
break
if params is None:
raise MissingParametersError(
f"Missing HarmonicAngle parameters for angles between "
f"atoms {atom_1_type}, {atom_2_type} and {atom_3_type}"
)
return params
def _extract_periodic_proper_params(self, atom_types):
"""Return parameters for a PeriodicTorsion (proper) between atom types."""
if len(atom_types) != 4:
raise ValueError(
f"PeriodicTorsion proper parameters can only "
f"be extracted for four atoms. Provided {len(atom_types)}"
)
periodic_torsion_force_gen = self.get_generator(
ff=self, gen_type=PeriodicTorsionGenerator
)
wildcard = self._atomClasses[""]
(
atom_1_type,
atom_2_type,
atom_3_type,
atom_4_type,
) = self.substitute_wildcards(atom_types, wildcard)
match = None
for index in periodic_torsion_force_gen.propersForAtomType[atom_2_type]:
tordef = periodic_torsion_force_gen.proper[index]
types1 = tordef.types1
types2 = tordef.types2
types3 = tordef.types3
types4 = tordef.types4
if (
atom_2_type in types2
and atom_3_type in types3
and atom_4_type in types4
and atom_1_type in types1
) or (
atom_2_type in types3
and atom_3_type in types2
and atom_4_type in types1
and atom_1_type in types4
):
has_wildcard = wildcard in (types1, types2, types3, types4)
if match is None or not has_wildcard:
match = tordef
if not has_wildcard:
break
if match is not None:
return self._get_periodic_torsion_params(match)
else:
raise MissingParametersError(
f"Missing PeriodicTorsion parameters for proper dihedrals between "
f"atoms {atom_types[0]}, {atom_types[1]}, {atom_types[2]} and {atom_types[3]}."
)
def _extract_periodic_improper_params(self, atom_types):
"""Return parameters for a PeriodicTorsion (improper) between atom types."""
if len(atom_types) != 4:
raise ValueError(
f"PeriodicTorsion improper parameters can only "
f"be extracted for four atoms. Provided {len(atom_types)}"
)
periodic_torsion_force_gen = self.get_generator(
ff=self, gen_type=PeriodicTorsionGenerator
)
match = self._match_impropers(atom_types, periodic_torsion_force_gen)
if match is not None:
return self._get_periodic_torsion_params(match)
else:
raise MissingParametersError(
f"Missing PeriodicTorsion parameters for improper dihedrals between "
f"atoms {atom_types[0]}, {atom_types[1]}, {atom_types[2]} and {atom_types[3]}."
)
def _extract_rb_proper_params(self, atom_types):
"""Return parameters for a RBTorsion (proper) between atom types."""
if len(atom_types) != 4:
raise ValueError(
f"RBTorsion proper parameters can only "
f"be extracted for four atoms. Provided {len(atom_types)}"
)
rb_torsion_force_gen = self.get_generator(
ff=self, gen_type=RBTorsionGenerator
)
wildcard = self._atomClasses[""]
(
atom_1_type,
atom_2_type,
atom_3_type,
atom_4_type,
) = self.substitute_wildcards(atom_types, wildcard)
match = None
for tordef in rb_torsion_force_gen.proper:
types1 = tordef.types1
types2 = tordef.types2
types3 = tordef.types3
types4 = tordef.types4
if (
atom_2_type in types2
and atom_3_type in types3
and atom_4_type in types4
and atom_1_type in types1
) or (
atom_2_type in types3
and atom_3_type in types2
and atom_4_type in types1
and atom_1_type in types4
):
has_wildcard = wildcard in (types1, types2, types3, types4)
if match is None or not has_wildcard:
match = tordef
if not has_wildcard:
break
if match is not None:
return self._get_rb_torsion_params(match)
else:
raise MissingParametersError(
f"Missing RBTorsion parameters for proper dihedrals between "
f"atoms {atom_types[0]}, {atom_types[1]}, {atom_types[2]} and {atom_types[3]}."
)
def _extract_rb_improper_params(self, atom_types):
"""Return parameters for a RBTorsion (improper) between atom types."""
if len(atom_types) != 4:
raise ValueError(
f"RBTorsion improper parameters can only "
f"be extracted for four atoms. Provided {len(atom_types)}"
)
rb_torsion_force_gen = self.get_generator(
ff=self, gen_type=RBTorsionGenerator
)
match = self._match_impropers(atom_types, rb_torsion_force_gen)
if match is not None:
return self._get_rb_torsion_params(match)
else:
raise MissingParametersError(
f"Missing RBTorsion parameters for improper dihedrals between "
f"atoms {atom_types[0]}, {atom_types[1]}, {atom_types[2]} and {atom_types[3]}."
)
[docs] def map_atom_classes_to_types(self, atom_classes_keys, strict=False):
"""Replace the atom_classes in the provided list by atom_types."""
atom_type_keys = []
for key in atom_classes_keys:
# When to do this substitution with wildcards?
substitution = self._atomClasses.get(key)
if not substitution:
raise ValueError(
f"Atom class {key} is missing from the Forcefield"
)
atom_type_keys.append(next(iter(substitution)))
return atom_type_keys
@staticmethod
def _match_impropers(atom_types, generator):
"""Find a match between an improper in the generator and atom types provided."""
wildcard = generator.ff._atomClasses[""]
(
atom_1_type,
atom_2_type,
atom_3_type,
atom_4_type,
) = Forcefield.substitute_wildcards(atom_types, wildcard)
match = None
for improper in generator.improper:
types1 = improper.types1
types2 = improper.types2
types3 = improper.types3
types4 = improper.types4
has_wildcard = wildcard in (types1, types2, types3, types4)
if match is not None and has_wildcard:
# Prefer specific definitions over ones with wildcards
continue
if atom_1_type in types1:
for (t2, t3, t4) in itertools.permutations(
((atom_2_type, 1), (atom_3_type, 2), (atom_4_type, 3))
):
if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
match = improper
break
return match
@staticmethod
def _get_rb_torsion_params(torsion):
return {
"c0": torsion.c[0],
"c1": torsion.c[1],
"c2": torsion.c[2],
"c3": torsion.c[3],
"c4": torsion.c[4],
"c5": torsion.c[5],
}
@staticmethod
def _get_periodic_torsion_params(torsion):
params = {"periodicity": [], "phase": [], "k": []}
for i in range(len(torsion.phase)):
params["periodicity"].append(torsion.periodicity[i])
params["phase"].append(torsion.phase[i])
params["k"].append(torsion.k[i])
return params
[docs] @staticmethod
def get_generator(ff, gen_type):
"""Return a specific force generator for this Forcefield.
Parameters
----------
ff: foyer.Forcefield
The forcefield to return generator for
gen_type: Type
The generator type to return
Returns
-------
instance of gen_type
The instance of `gen_type` for this Forcefield
Raises
------
MissingForceError
If the Forcefield doesn't have a generator of type `gen_type`
"""
generator = list(
filter(
lambda x: isinstance(x, gen_type),
ff.getGenerators(),
)
)
if generator:
return generator.pop()
else:
raise MissingForceError(
f"{force_for[gen_type]} is missing from the Forcefield"
)
[docs] @staticmethod
def substitute_wildcards(atom_types, wildcard):
"""Return possible wildcard options."""
return tuple(
atom_type or next(iter(wildcard)) for atom_type in atom_types
)
pmd.Structure.write_foyer = write_foyer