6. Creating Force Fields

OpenMM uses a simple XML file format to describe force fields. It includes many common force fields, but you can also create your own. A force field can use all the standard OpenMM force classes, as well as the very flexible custom force classes. You can even extend the ForceField class to add support for completely new forces, such as ones defined in plugins. This makes it a powerful tool for force field development.

6.1. Basic Concepts

Let’s start by considering how OpenMM defines a force field. There are a small number of basic concepts to understand.

6.1.1. Atom Types and Atom Classes

Force field parameters are assigned to atoms based on their “atom types”. Atom types should be the most specific identification of an atom that will ever be needed. Two atoms should have the same type only if the force field will always treat them identically in every way.

Multiple atom types can be grouped together into “atom classes”. In general, two types should be in the same class if the force field usually (but not necessarily always) treats them identically. For example, the \(\alpha\)-carbon of an alanine residue will probably have a different atom type than the \(\alpha\)-carbon of a leucine residue, but both of them will probably have the same atom class.

All force field parameters can be specified either by atom type or atom class. Classes exist as a convenience to make force field definitions more compact. If necessary, you could define everything in terms of atom types, but when many types all share the same parameters, it is convenient to only have to specify them once.

6.1.2. Residue Templates

Types are assigned to atoms by matching residues to templates. A template specifies a list of atoms, the type of each one, and the bonds between them. For each residue in the PDB file, the force field searches its list of templates for one that has an identical set of atoms with identical bonds between them. When matching templates, neither the order of the atoms nor their names matter; it only cares about their elements and the set of bonds between them. (The PDB file reader does care about names, of course, since it needs to figure out which atom each line of the file corresponds to.)

6.1.3. Forces

Once a force field has defined its atom types and residue templates, it must define its force field parameters. This generally involves one block of XML for each Force object that will be added to the System. The details are different for each Force, but it generally consists of a set of rules for adding interactions based on bonds and atom types or classes. For example, when adding a HarmonicBondForce, the force field will loop over every pair of bonded atoms, check their types and classes, and see if they match any of its rules. If so, it will call addBond() on the HarmonicBondForce. If none of them match, it simply ignores that pair and continues.

6.2. Writing the XML File

The root element of the XML file must be a <ForceField> tag:

<ForceField>
...
</ForceField>

The <ForceField> tag contains the following children:

  • An <AtomTypes> tag containing the atom type definitions

  • A <Residues> tag containing the residue template definitions

  • Zero or more tags defining specific forces

The order of these tags does not matter. They are described in detail below.

6.2.1. <AtomTypes>

The atom type definitions look like this:

<AtomTypes>
 <Type name="0" class="N" element="N" mass="14.00672"/>
 <Type name="1" class="H" element="H" mass="1.007947"/>
 <Type name="2" class="CT" element="C" mass="12.01078"/>
 ...
</AtomTypes>

There is one <Type> tag for each atom type. It specifies the name of the type, the name of the class it belongs to, the symbol for its element, and its mass in amu. The names are arbitrary strings: they need not be numbers, as in this example. The only requirement is that all types have unique names. The classes are also arbitrary strings, and in general will not be unique. Two types belong to the same class if they list the same value for the class attribute.

6.2.2. <Residues>

The residue template definitions look like this:

<Residues>
 <Residue name="ACE">
  <Atom name="HH31" type="710"/>
  <Atom name="CH3" type="711"/>
  <Atom name="HH32" type="710"/>
  <Atom name="HH33" type="710"/>
  <Atom name="C" type="712"/>
  <Atom name="O" type="713"/>
  <Bond atomName1="HH31" atomName2="CH3"/>
  <Bond atomName1="CH3" atomName2="HH32"/>
  <Bond atomName1="CH3" atomName2="HH33"/>
  <Bond atomName1="CH3" atomName2="C"/>
  <Bond atomName1="C" atomName2="O"/>
  <ExternalBond atomName="C"/>
 </Residue>
 <Residue name="ALA">
  ...
 </Residue>
 ...
</Residues>

There is one <Residue> tag for each residue template. That in turn contains the following tags:

  • An <Atom> tag for each atom in the residue. This specifies the name of the atom and its atom type.

  • A <Bond> tag for each pair of atoms that are bonded to each other. The atomName1 and atomName2 attributes are the names of the two bonded atoms. (Some older force fields use the alternate tags to and from to specify the atoms by index instead of name. This is still supported for backward compatibility, but specifying atoms by name is recommended, since it makes the residue definition much easier to understand.)

  • An <ExternalBond> tag for each atom that will be bonded to an atom of a different residue. atomName is the name of the atom. (Alternatively, the deprecated from tag may indicate the atom by index instead of name.)

The <Residue> tag may also contain <VirtualSite> tags, as in the following example:

<Residue name="HOH">
 <Atom name="O" type="tip4pew-O"/>
 <Atom name="H1" type="tip4pew-H"/>
 <Atom name="H2" type="tip4pew-H"/>
 <Atom name="M" type="tip4pew-M"/>
 <VirtualSite type="average3" siteName="M" atomName1="O" atomName2="H1" atomName3="H2"
     weight1="0.786646558" weight2="0.106676721" weight3="0.106676721"/>
 <Bond atomName1="O" atomName2="H1"/>
 <Bond atomName1="O" atomName2="H2"/>
</Residue>

Each <VirtualSite> tag indicates an atom in the residue that should be represented with a virtual site. The type attribute may equal "average2", "average3", "outOfPlane", or "localCoords", which correspond to the TwoParticleAverageSite, ThreeParticleAverageSite, OutOfPlaneSite, and LocalCoordinatesSite classes respectively. The siteName attribute gives the name of the atom to represent with a virtual site. The atoms it is calculated based on are specified by atomName1, atomName2, etc. (Some old force fields use the deprecated tags index, atom1, atom2, etc. to refer to them by index instead of name.)

The remaining attributes are specific to the virtual site class, and specify the parameters for calculating the site position. For a TwoParticleAverageSite, they are weight1 and weight2. For a ThreeParticleAverageSite, they are weight1, weight2, and weight3. For an OutOfPlaneSite, they are weight12, weight13, and weightCross. For a LocalCoordinatesSite, they are p1, p2, and p3 (giving the x, y, and z coordinates of the site position in the local coordinate system), and wo1, wx1, wy1, wo2, wx2, wy2, … (giving the weights for computing the origin, x axis, and y axis).

6.2.3. <Patches>

A “patch” is a set of rules for modifying a residue template (or possibly multiple templates at once). For example a terminal amino acid is slightly different from one in the middle of a chain. A force field could of course define multiple templates for each amino acid (standard, N-terminal, C-terminal, and monomer), but since the modifications are the same for nearly all amino acids, it is simpler to include only the “standard” templates, along with a set of patches for modifying terminal residues.

Here is an example of a patch definition:

<Patch name="NTER">
 <RemoveAtom name="H"/>
 <RemoveBond atomName1="N" atomName2="H"/>
 <AddAtom name="H1" type="H"/>
 <AddAtom name="H2" type="H"/>
 <AddAtom name="H3" type="H"/>
 <AddBond atomName1="N" atomName2="H1"/>
 <AddBond atomName1="N" atomName2="H2"/>
 <AddBond atomName1="N" atomName2="H3"/>
 <RemoveExternalBond atomName="N"/>
 <ChangeAtom name="N" type="N3"/>
</Patch>

There is one <Patch> tag for each patch definition. That in turn may contain any of the following tags:

  • An <AddAtom> tag indicates that an atom should be added to the template. It specifies the name of the atom and its atom type.

  • A <ChangeAtom> tag indicates that the type of an atom already present in the template should be altered. It specifies the name of the atom and its new atom type.

  • A <RemoveAtom> tag indicates that an atom should be removed from the template. It specifies the name of the atom to remove.

  • An <AddBond> tag indicates that a bond should be added to the template. It specifies the names of the two bonded atoms.

  • A <RemoveBond> tag indicates that a bond already present in the template should be removed. It specifies the names of the two bonded atoms.

  • An <AddExternalBond> tag indicates that a new external bond should be added to the template. It specifies the name of the bonded atom.

  • A <RemoveExternalBond> tag indicates that an external bond aleady present in the template should be removed. It specifies the name of the bonded atom.

In addition to defining the patches, you also must identify which residue templates each patch can be applied to. This can be done in two ways. The more common one is to have each template identify the patches that can be applied to it. This is done with an <AllowPatch> tag:

<Residue name="ALA">
 <AllowPatch name="CTER"/>
 <AllowPatch name="NTER"/>
 ...
</Residue>

Alternatively, the patch can indicate which residues it may be applied to. This is done with an <ApplyToResidue> tag:

<Patch name="NTER">
 <ApplyToResidue name="ALA"/>
 <ApplyToResidue name="ARG"/>
 ...
</Patch>

A patch can alter multiple templates at once. This is useful for creating bonds between molecules, and allows the atom types in one residue to depend on the identity of the other residue it is bonded to. To create a multi-residue patch, added a residues attribute to the <Patch> tag specifying how many residues that patch covers. Then whenever you refer to an atom, prefix its name with the index of the residue it belongs to:

<Patch name="Disulfide" residues="2">
  <RemoveAtom name="1:HG"/>
  <RemoveAtom name="2:HG"/>
  <AddBond atomName1="1:SG" atomName2="2:SG"/>
  <ApplyToResidue name="1:CYS"/>
  <ApplyToResidue name="2:CYS"/>
</Patch>

In this example, the patch modifies two residues of the same type, but that need not always be true. Each <ApplyToResidue> tag therefore indicates which one of the residue templates it modifies may be of the specified type. Similarly, if a residue template includes an <AcceptPatch> tag for a multi-residue patch, it must specify the name of the patch, followed by the index of the residue within that patch:

<AllowPatch name="Disulfide:1"/>

6.2.4. Missing residue templates

Caution

These features are experimental, and their API is subject to change.

You can use the getUnmatchedResidues() method to get a list of residues in the provided topology object that do not currently have a matching residue template defined in the ForceField.

pdb = PDBFile('input.pdb')
forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
unmatched_residues = forcefield.getUnmatchedResidues(topology)

This is useful for identifying issues with prepared systems, debugging issues with residue template definitions, or identifying which additional residues need to be parameterized.

As a convenience for parameterizing new residues, you can also get a list of residues and empty residue templates using generateTemplatesForUnmatchedResidues()

pdb = PDBFile('input.pdb')
forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
[templates, residues] = forcefield.generateTemplatesForUnmatchedResidues(topology)
# Se the atom types
for template in templates:
    for atom in template.atoms:
        atom.type = ... # set the atom types here
    # Register the template with the forcefield.
    forcefield.registerResidueTemplate(template)

If you find that templates seem to be incorrectly matched, another useful function getMatchingTemplates() can help you identify which templates are being matched:

pdb = PDBFile('input.pdb')
forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
templates = forcefield.getMatchingTemplates(topology)
for (residue, template) in zip(pdb.topology.residues(), templates):
    print("Residue %d %s matched template %s" % (residue.id, residue.name, template.name))

6.2.5. <HarmonicBondForce>

To add a HarmonicBondForce to the System, include a tag that looks like this:

<HarmonicBondForce>
 <Bond class1="C" class2="C" length="0.1525" k="259408.0"/>
 <Bond class1="C" class2="CA" length="0.1409" k="392459.2"/>
 <Bond class1="C" class2="CB" length="0.1419" k="374049.6"/>
 ...
</HarmonicBondForce>

Every <Bond> tag defines a rule for creating harmonic bond interactions between atoms. Each tag may identify the atoms either by type (using the attributes type1 and type2) or by class (using the attributes class1 and class2). For every pair of bonded atoms, the force field searches for a rule whose atom types or atom classes match the two atoms. If it finds one, it calls addBond() on the HarmonicBondForce with the specified parameters. Otherwise, it ignores that pair and continues. length is the equilibrium bond length in nm, and k is the spring constant in kJ/mol/nm2.

6.2.6. <HarmonicAngleForce>

To add a HarmonicAngleForce to the System, include a tag that looks like this:

<HarmonicAngleForce>
 <Angle class1="C" class2="C" class3="O" angle="2.094" k="669.44"/>
 <Angle class1="C" class2="C" class3="OH" angle="2.094" k="669.44"/>
 <Angle class1="CA" class2="C" class3="CA" angle="2.094" k="527.184"/>
 ...
</HarmonicAngleForce>

Every <Angle> tag defines a rule for creating harmonic angle interactions between triplets of atoms. Each tag may identify the atoms either by type (using the attributes type1, type2, …) or by class (using the attributes class1, class2, …). The force field identifies every set of three atoms in the system where the first is bonded to the second, and the second to the third. For each one, it searches for a rule whose atom types or atom classes match the three atoms. If it finds one, it calls addAngle() on the HarmonicAngleForce with the specified parameters. Otherwise, it ignores that set and continues. angle is the equilibrium angle in radians, and k is the spring constant in kJ/mol/radian2.

6.2.7. <PeriodicTorsionForce>

To add a PeriodicTorsionForce to the System, include a tag that looks like this:

<PeriodicTorsionForce>
 <Proper class1="HC" class2="CT" class3="CT" class4="CT" periodicity1="3" phase1="0.0"
     k1="0.66944"/>
 <Proper class1="HC" class2="CT" class3="CT" class4="HC" periodicity1="3" phase1="0.0"
     k1="0.6276"/>
 ...
 <Improper class1="N" class2="C" class3="CT" class4="O" periodicity1="2"
     phase1="3.14159265359" k1="4.6024"/>
 <Improper class1="N" class2="C" class3="CT" class4="H" periodicity1="2"
     phase1="3.14159265359" k1="4.6024"/>
 ...
</PeriodicTorsionForce>

Every child tag defines a rule for creating periodic torsion interactions between sets of four atoms. Each tag may identify the atoms either by type (using the attributes type1, type2, …) or by class (using the attributes class1, class2, …).

The force field recognizes two different types of torsions: proper and improper. A proper torsion involves four atoms that are bonded in sequence: 1 to 2, 2 to 3, and 3 to 4. An improper torsion involves a central atom and three others that are bonded to it: atoms 2, 3, and 4 are all bonded to atom 1. The force field begins by identifying every set of atoms in the system of each of these types. For each one, it searches for a rule whose atom types or atom classes match the four atoms. If it finds one, it calls addTorsion() on the PeriodicTorsionForce with the specified parameters. Otherwise, it ignores that set and continues. periodicity1 is the periodicity of the torsion, phase1 is the phase offset in radians, and k1 is the force constant in kJ/mol.

Each torsion definition can specify multiple periodic torsion terms to add to its atoms. To add a second one, just add three more attributes: periodicity2, phase2, and k2. You can have as many terms as you want. Here is an example of a rule that adds three torsion terms to its atoms:

<Proper class1="CT" class2="CT" class3="CT" class4="CT"
    periodicity1="3" phase1="0.0" k1="0.75312"
    periodicity2="2" phase2="3.14159265359" k2="1.046"
    periodicity3="1" phase3="3.14159265359" k3="0.8368"/>

You can also use wildcards when defining torsions. To do this, simply leave the type or class name for an atom empty. That will cause it to match any atom. For example, the following definition will match any sequence of atoms where the second atom has class OS and the third has class P:

<Proper class1="" class2="OS" class3="P" class4="" periodicity1="3" phase1="0.0" k1="1.046"/>

The <PeriodicTorsionForce> tag also supports an optional ordering attribute to provide better compatibility with the way impropers are assigned in different simulation packages:

  • ordering="default" specifies the default behavior if the attribute is omitted.

  • ordering="amber" produces behavior that replicates the behavior of AmberTools LEaP

  • ordering="charmm" produces behavior more consistent with CHARMM

  • ordering="smirnoff" allows multiple impropers to be added using exact matching to replicate the beheavior of SMIRNOFF impropers

Different <PeriodicTorsionForce> tags can specify different ordering values to be used for the sub-elements appearing within their tags.

6.2.8. <RBTorsionForce>

To add an RBTorsionForce to the System, include a tag that looks like this:

<RBTorsionForce>
 <Proper class1="CT" class2="CT" class3="OS" class4="CT" c0="2.439272" c1="4.807416"
     c2="-0.8368" c3="-6.409888" c4="0" c5="0" />
 <Proper class1="C" class2="N" class3="CT" class4="C" c0="10.46" c1="-3.34720"
     c2="-7.1128" c3="0" c4="0" c5="0" />
 ...
 <Improper class1="N" class2="C" class3="CT" class4="O" c0="0.8368" c1="0"
     c2="-2.76144" c3="0" c4="3.3472" c5="0" />
 <Improper class1="N" class2="C" class3="CT" class4="H" c0="29.288" c1="-8.368"
     c2="-20.92" c3="0" c4="0" c5="0" />
 ...
</RBTorsionForce>

Every child tag defines a rule for creating Ryckaert-Bellemans torsion interactions between sets of four atoms. Each tag may identify the atoms either by type (using the attributes type1, type2, …) or by class (using the attributes class1, class2, …).

The force field recognizes two different types of torsions: proper and improper. A proper torsion involves four atoms that are bonded in sequence: 1 to 2, 2 to 3, and 3 to 4. An improper torsion involves a central atom and three others that are bonded to it: atoms 2, 3, and 4 are all bonded to atom 1. The force field begins by identifying every set of atoms in the system of each of these types. For each one, it searches for a rule whose atom types or atom classes match the four atoms. If it finds one, it calls addTorsion() on the RBTorsionForce with the specified parameters. Otherwise, it ignores that set and continues. The attributes c0 through c5 are the coefficients of the terms in the Ryckaert-Bellemans force expression.

You can also use wildcards when defining torsions. To do this, simply leave the type or class name for an atom empty. That will cause it to match any atom. For example, the following definition will match any sequence of atoms where the second atom has class OS and the third has class P:

<Proper class1="" class2="OS" class3="P" class4="" c0="2.439272" c1="4.807416"
    c2="-0.8368" c3="-6.409888" c4="0" c5="0" />

6.2.9. <CMAPTorsionForce>

To add a CMAPTorsionForce to the System, include a tag that looks like this:

<CMAPTorsionForce>
 <Map>
  0.0 0.809 0.951 0.309
  -0.587 -1.0 -0.587 0.309
  0.951 0.809 0.0 -0.809
  -0.951 -0.309 0.587 1.0
 </Map>
 <Torsion map="0" class1="CT" class2="CT" class3="C" class4="N" class5="CT"/>
 <Torsion map="0" class1="N" class2="CT" class3="C" class4="N" class5="CT"/>
 ...
</CMAPTorsionForce>

Each <Map> tag defines an energy correction map. Its content is the list of energy values in kJ/mole, listed in the correct order for CMAPTorsionForce’s addMap() method and separated by white space. See the API documentation for details. The size of the map is determined from the number of energy values.

Each <Torsion> tag defines a rule for creating CMAP torsion interactions between sets of five atoms. The tag may identify the atoms either by type (using the attributes type1, type2, …) or by class (using the attributes class1, class2, …). The force field identifies every set of five atoms that are bonded in sequence: 1 to 2, 2 to 3, 3 to 4, and 4 to 5. For each one, it searches for a rule whose atom types or atom classes match the five atoms. If it finds one, it calls addTorsion() on the CMAPTorsionForce with the specified parameters. Otherwise, it ignores that set and continues. The first torsion is defined by the sequence of atoms 1-2-3-4, and the second one by atoms 2-3-4-5. map is the index of the map to use, starting from 0, in the order they are listed in the file.

You can also use wildcards when defining torsions. To do this, simply leave the type or class name for an atom empty. That will cause it to match any atom. For example, the following definition will match any sequence of five atoms where the middle three have classes CT, C, and N respectively:

<Torsion map="0" class1="" class2="CT" class3="C" class4="N" class5=""/>

6.2.10. <NonbondedForce>

To add a NonbondedForce to the System, include a tag that looks like this:

<NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
 <Atom type="0" charge="-0.4157" sigma="0.32499" epsilon="0.71128"/>
 <Atom type="1" charge="0.2719" sigma="0.10690" epsilon="0.06568"/>
 <Atom type="2" charge="0.0337" sigma="0.33996" epsilon="0.45772"/>
 ...
</NonbondedForce>

The <NonbondedForce> tag has two attributes coulomb14scale and lj14scale that specify the scale factors between pairs of atoms separated by three bonds. After setting the nonbonded parameters for all atoms, the force field calls createExceptionsFromBonds() on the NonbondedForce, passing in these scale factors as arguments.

Each <Atom> tag specifies the nonbonded parameters for one atom type (specified with the type attribute) or atom class (specified with the class attribute). It is fine to mix these two methods, having some tags specify a type and others specify a class. However you do it, you must make sure that a unique set of parameters is defined for every atom type. charge is measured in units of the proton charge, sigma is in nm, and epsilon is in kJ/mole.

6.2.11. <GBSAOBCForce>

To add a GBSAOBCForce to the System, include a tag that looks like this:

<GBSAOBCForce>
 <Atom type="0" charge="-0.4157" radius="0.1706" scale="0.79"/>
 <Atom type="1" charge="0.2719" radius="0.115" scale="0.85"/>
 <Atom type="2" charge="0.0337" radius="0.19" scale="0.72"/>
 ...
</GBSAOBCForce>

Each <Atom> tag specifies the OBC parameters for one atom type (specified with the type attribute) or atom class (specified with the class attribute). It is fine to mix these two methods, having some tags specify a type and others specify a class. However you do it, you must make sure that a unique set of parameters is defined for every atom type. charge is measured in units of the proton charge, radius is the GBSA radius in nm, and scale is the OBC scaling factor.

6.2.12. <CustomBondForce>

To add a CustomBondForce to the System, include a tag that looks like this:

<CustomBondForce energy="scale*k*(r-r0)^2">
 <GlobalParameter name="scale" defaultValue="0.5"/>
 <PerBondParameter name="k"/>
 <PerBondParameter name="r0"/>
 <Bond class1="OW" class2="HW" r0="0.09572" k="462750.4"/>
 <Bond class1="HW" class2="HW" r0="0.15136" k="462750.4"/>
 <Bond class1="C" class2="C" r0="0.1525" k="259408.0"/>
 ...
</CustomBondForce>

The energy expression for the CustomBondForce is specified by the energy attribute. This is a mathematical expression that gives the energy of each bond as a function of its length r. It also may depend on an arbitrary list of global or per-bond parameters. Use a <GlobalParameter> tag to define a global parameter, and a <PerBondParameter> tag to define a per-bond parameter.

Every <Bond> tag defines a rule for creating custom bond interactions between atoms. Each tag may identify the atoms either by type (using the attributes type1 and type2) or by class (using the attributes class1 and class2). For every pair of bonded atoms, the force field searches for a rule whose atom types or atom classes match the two atoms. If it finds one, it calls addBond() on the CustomBondForce. Otherwise, it ignores that pair and continues. The remaining attributes are the values to use for the per-bond parameters. All per-bond parameters must be specified for every <Bond> tag, and the attribute name must match the name of the parameter. For instance, if there is a per-bond parameter with the name “k”, then every <Bond> tag must include an attribute called k.

6.2.13. <CustomAngleForce>

To add a CustomAngleForce to the System, include a tag that looks like this:

<CustomAngleForce energy="scale*k*(theta-theta0)^2">
 <GlobalParameter name="scale" defaultValue="0.5"/>
 <PerAngleParameter name="k"/>
 <PerAngleParameter name=" theta0"/>
 <Angle class1="HW" class2="OW" class3="HW" theta0="1.824218" k="836.8"/>
 <Angle class1="HW" class2="HW" class3="OW" theta0="2.229483" k="0.0"/>
 <Angle class1="C" class2="C" class3="O" theta0="2.094395" k="669.44"/>
 ...
</CustomAngleForce>

The energy expression for the CustomAngleForce is specified by the energy attribute. This is a mathematical expression that gives the energy of each angle as a function of the angle theta. It also may depend on an arbitrary list of global or per-angle parameters. Use a <GlobalParameter> tag to define a global parameter, and a <PerAngleParameter> tag to define a per-angle parameter.

Every <Angle> tag defines a rule for creating custom angle interactions between triplets of atoms. Each tag may identify the atoms either by type (using the attributes type1, type2, …) or by class (using the attributes class1, class2, …). The force field identifies every set of three atoms in the system where the first is bonded to the second, and the second to the third. For each one, it searches for a rule whose atom types or atom classes match the three atoms. If it finds one, it calls addAngle() on the CustomAngleForce. Otherwise, it ignores that set and continues. The remaining attributes are the values to use for the per-angle parameters. All per-angle parameters must be specified for every <Angle> tag, and the attribute name must match the name of the parameter. For instance, if there is a per-angle parameter with the name “k”, then every <Angle> tag must include an attribute called k.

6.2.14. <CustomTorsionForce>

To add a CustomTorsionForce to the System, include a tag that looks like this:

<CustomTorsionForce energy="scale*k*(1+cos(per*theta-phase))">
 <GlobalParameter name="scale" defaultValue="1"/>
 <PerTorsionParameter name="k"/>
 <PerTorsionParameter name="per"/>
 <PerTorsionParameter name="phase"/>
 <Proper class1="HC" class2="CT" class3="CT" class4="CT" per="3" phase="0.0" k="0.66944"/>
 <Proper class1="HC" class2="CT" class3="CT" class4="HC" per="3" phase="0.0" k="0.6276"/>
 ...
 <Improper class1="N" class2="C" class3="CT" class4="O" per="2" phase="3.14159265359"
     k="4.6024"/>
 <Improper class1="N" class2="C" class3="CT" class4="H" per="2" phase="3.14159265359"
     k="4.6024"/>
 ...
</CustomTorsionForce>

The energy expression for the CustomTorsionForce is specified by the energy attribute. This is a mathematical expression that gives the energy of each torsion as a function of the angle theta. It also may depend on an arbitrary list of global or per-torsion parameters. Use a <GlobalParameter> tag to define a global parameter, and a <PerTorsionParameter> tag to define a per-torsion parameter.

Every child tag defines a rule for creating custom torsion interactions between sets of four atoms. Each tag may identify the atoms either by type (using the attributes type1, type2, …) or by class (using the attributes class1, class2, …).

The force field recognizes two different types of torsions: proper and improper. A proper torsion involves four atoms that are bonded in sequence: 1 to 2, 2 to 3, and 3 to 4. An improper torsion involves a central atom and three others that are bonded to it: atoms 2, 3, and 4 are all bonded to atom 1. The force field begins by identifying every set of atoms in the system of each of these types. For each one, it searches for a rule whose atom types or atom classes match the four atoms. If it finds one, it calls addTorsion() on the CustomTorsionForce with the specified parameters. Otherwise, it ignores that set and continues. The remaining attributes are the values to use for the per- torsion parameters. Every <Torsion> tag must include one attribute for every per-torsion parameter, and the attribute name must match the name of the parameter.

You can also use wildcards when defining torsions. To do this, simply leave the type or class name for an atom empty. That will cause it to match any atom. For example, the following definition will match any sequence of atoms where the second atom has class OS and the third has class P:

<Proper class1="" class2="OS" class3="P" class4="" per="3" phase="0.0" k="0.66944"/>

6.2.15. <CustomNonbondedForce>

To add a CustomNonbondedForce to the System, include a tag that looks like this:

<CustomNonbondedForce energy="scale*epsilon1*epsilon2*((sigma1+sigma2)/r)^12" bondCutoff="3">
 <GlobalParameter name="scale" defaultValue="1"/>
 <PerParticleParameter name="sigma"/>
 <PerParticleParameter name="epsilon"/>
 <Atom type="0" sigma="0.3249" epsilon="0.7112"/>
 <Atom type="1" sigma="0.1069" epsilon="0.0656"/>
 <Atom type="2" sigma="0.3399" epsilon="0.4577"/>
 ...
</CustomNonbondedForce>

The energy expression for the CustomNonbondedForce is specified by the energy attribute. This is a mathematical expression that gives the energy of each pairwise interaction as a function of the distance r. It also may depend on an arbitrary list of global or per-particle parameters. Use a <GlobalParameter> tag to define a global parameter, and a <PerParticleParameter> tag to define a per-particle parameter.

Exclusions are created automatically based on the bondCutoff attribute. After setting the nonbonded parameters for all atoms, the force field calls createExclusionsFromBonds() on the CustomNonbondedForce, passing in this value as its argument. To avoid creating exclusions, set bondCutoff to 0.

Each <Atom> tag specifies the parameters for one atom type (specified with the type attribute) or atom class (specified with the class attribute). It is fine to mix these two methods, having some tags specify a type and others specify a class. However you do it, you must make sure that a unique set of parameters is defined for every atom type. The remaining attributes are the values to use for the per-atom parameters. All per-atom parameters must be specified for every <Atom> tag, and the attribute name must match the name of the parameter. For instance, if there is a per-atom parameter with the name “radius”, then every <Atom> tag must include an attribute called radius.

CustomNonbondedForce also allows you to define tabulated functions. See Section 6.2.20 for details.

6.2.16. <CustomGBForce>

To add a CustomGBForce to the System, include a tag that looks like this:

<CustomGBForce>
 <GlobalParameter name="solventDielectric" defaultValue="78.3"/>
 <GlobalParameter name="soluteDielectric" defaultValue="1"/>
 <PerParticleParameter name="charge"/>
 <PerParticleParameter name="radius"/>
 <PerParticleParameter name="scale"/>
 <ComputedValue name="I" type="ParticlePairNoExclusions">
    step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);
    U=r+sr2; C=2*(1/or1-1/L)*step(sr2-r-or1); L=max(or1, D); D=abs(r-sr2); sr2 =
    scale2*or2; or1 = radius1-0.009; or2 = radius2-0.009
 </ComputedValue>
 <ComputedValue name="B" type="SingleParticle">
  1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius); psi=I*or; or=radius-0.009
 </ComputedValue>
 <EnergyTerm type="SingleParticle">
  28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935456*
          (1/soluteDielectric-1/solventDielectric)*charge^2/B
 </EnergyTerm>
 <EnergyTerm type="ParticlePair">
  -138.935456*(1/soluteDielectric-1/solventDielectric)*charge1*charge2/f;
          f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))
 </EnergyTerm>
 <Atom type="0" charge="-0.4157" radius="0.1706" scale="0.79"/>
 <Atom type="1" charge="0.2719" radius="0.115" scale="0.85"/>
 <Atom type="2" charge="0.0337" radius="0.19" scale="0.72"/>
 ...
</CustomGBForce>

The above (rather complicated) example defines a generalized Born model that is equivalent to GBSAOBCForce. The definition consists of a set of computed values (defined by <ComputedValue> tags) and energy terms (defined by <EnergyTerm> tags), each of which is evaluated according to a mathematical expression. See the API documentation for details.

The expressions may depend on an arbitrary list of global or per-atom parameters. Use a <GlobalParameter> tag to define a global parameter, and a <PerAtomParameter> tag to define a per-atom parameter.

Each <Atom> tag specifies the parameters for one atom type (specified with the type attribute) or atom class (specified with the class attribute). It is fine to mix these two methods, having some tags specify a type and others specify a class. However you do it, you must make sure that a unique set of parameters is defined for every atom type. The remaining attributes are the values to use for the per-atom parameters. All per-atom parameters must be specified for every <Atom> tag, and the attribute name must match the name of the parameter. For instance, if there is a per-atom parameter with the name “radius”, then every <Atom> tag must include an attribute called radius.

CustomGBForce also allows you to define tabulated functions. See Section 6.2.20 for details.

6.2.17. <CustomHbondForce>

To add a CustomHbondForce to the System, include a tag that looks like this:

<CustomHbondForce particlesPerDonor="3" particlesPerAcceptor="2" bondCutoff="2"
    energy="scale*k*(distance(a1,d1)-r0)^2*(angle(a1,d1,d2)-theta0)^2">
 <GlobalParameter name="scale" defaultValue="1"/>
 <PerDonorParameter name="theta0"/>
 <PerAcceptorParameter name="k"/>
 <PerAcceptorParameter name="r0"/>
 <Donor class1="H" class2="N" class3="C" theta0="2.1"/>
 <Acceptor class1="O" class2="C" k="115.0" r0="0.2"/>
 ...
</CustomHbondForce>

The energy expression for the CustomHbondForce is specified by the energy attribute. This is a mathematical expression that gives the energy of each donor-acceptor interaction as a function of various particle coordinates, distances, and angles. See the API documentation for details. particlesPerDonor specifies the number of particles that make up a donor group, and particlesPerAcceptor specifies the number of particles that make up an acceptor group.

The expression may depend on an arbitrary list of global, per-donor, or per-acceptor parameters. Use a <GlobalParameter> tag to define a global parameter, a <PerDonorParameter> tag to define a per-donor parameter, and a <PerAcceptorParameter> tag to define a per-acceptor parameter.

Exclusions are created automatically based on the bondCutoff attribute. If any atom of a donor is within the specified distance (measured in bonds) of any atom of an acceptor, an exclusion is added to prevent them from interacting with each other. If a donor and an acceptor share any atom in common, that is a bond distance of 0, so they are always excluded.

Every <Donor> or <Acceptor> tag defines a rule for creating donor or acceptor groups. The number of atoms specified in each one must match the value of particlesPerDonor or particlesPerAcceptor specified in the parent tag. Each tag may identify the atoms either by type (using the attributes type1, type2, …) or by class (using the attributes class1, class2, …). The force field considers every atom in the system (if the number of atoms is 1), every pair of bonded atoms (if the number of atoms is 2), or every set of three atoms where the first is bonded to the second and the second to the third (if the number of atoms is 3). For each one, it searches for a rule whose atom types or atom classes match the atoms. If it finds one, it calls addDonor() or addAcceptor() on the CustomHbondForce. Otherwise, it ignores that set and continues. The remaining attributes are the values to use for the per-donor and per-acceptor parameters. All parameters must be specified for every tag, and the attribute name must match the name of the parameter. For instance, if there is a per-donor parameter with the name “k”, then every <Donor> tag must include an attribute called k.

CustomHbondForce also allows you to define tabulated functions. See Section 6.2.20 for details.

6.2.18. <CustomManyParticleForce>

To add a CustomManyParticleForce to the System, include a tag that looks like this:

<CustomManyParticleForce particlesPerSet="3" permutationMode="UniqueCentralParticle"
    bondCutoff="3" energy="scale*(distance(p1,p2)-r1)*(distance(p1,p3)-r1)">
 <GlobalParameter name="scale" defaultValue="1"/>
 <PerParticleParameter name="r"/>
 <TypeFilter index="0" types="1,2"/>
 <Atom type="0" r="0.31" filterType="0"/>
 <Atom type="1" r="0.25" filterType="0"/>
 <Atom type="2" r="0.33" filterType="1"/>
 ...
</CustomManyParticleForce>

The energy expression for the CustomManyParticleForce is specified by the energy attribute. This is a mathematical expression that gives the energy of each interaction as a function of various particle coordinates, distances, and angles. See the API documentation for details. particlesPerSet specifies the number of particles involved in the interaction and permutationMode specifies the permutation mode.

The expression may depend on an arbitrary list of global or per-atom parameters. Use a <GlobalParameter> tag to define a global parameter, and a <PerAtomParameter> tag to define a per-atom parameter.

Exclusions are created automatically based on the bondCutoff attribute. After setting the nonbonded parameters for all atoms, the force field calls createExclusionsFromBonds() on the CustomManyParticleForce, passing in this value as its argument. To avoid creating exclusions, set bondCutoff to 0.

Type filters may be specified with a <TypeFilter> tag. The index attribute specifies the index of the particle to apply the filter to, and types is a comma separated list of allowed types.

Each <Atom> tag specifies the parameters for one atom type (specified with the type attribute) or atom class (specified with the class attribute). It is fine to mix these two methods, having some tags specify a type and others specify a class. However you do it, you must make sure that a unique set of parameters is defined for every atom type. In addition, each <Atom> tag must include the filterType attribute, which specifies the atom type for use in type filters. The remaining attributes are the values to use for the per-atom parameters. All per-atom parameters must be specified for every <Atom> tag, and the attribute name must match the name of the parameter. For instance, if there is a per-atom parameter with the name “radius”, then every <Atom> tag must include an attribute called radius.

CustomManyParticleForce also allows you to define tabulated functions. See Section 6.2.20 for details.

6.2.19. Writing Custom Expressions

The custom forces described in this chapter involve user defined algebraic expressions. These expressions are specified as character strings, and may involve a variety of standard operators and mathematical functions.

The following operators are supported: + (add), - (subtract), * (multiply), / (divide), and ^ (power). Parentheses “(“and “)” may be used for grouping.

The following standard functions are supported: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select. step(x) = 0 if x < 0, 1 otherwise. delta(x) = 1 if x is 0, 0 otherwise. select(x,y,z) = z if x = 0, y otherwise. Some custom forces allow additional functions to be defined from tabulated values.

Numbers may be given in either decimal or exponential form. All of the following are valid numbers: 5, -3.1, 1e6, and 3.12e-2.

The variables that may appear in expressions are specified in the API documentation for each force class. In addition, an expression may be followed by definitions for intermediate values that appear in the expression. A semicolon “;” is used as a delimiter between value definitions. For example, the expression

a^2+a*b+b^2; a=a1+a2; b=b1+b2

is exactly equivalent to

(a1+a2)^2+(a1+a2)*(b1+b2)+(b1+b2)^2

The definition of an intermediate value may itself involve other intermediate values. All uses of a value must appear before that value’s definition.

6.2.20. Tabulated Functions

Some forces, such as CustomNonbondedForce and CustomGBForce, allow you to define tabulated functions. To define a function, include a <Function> tag inside the <CustomNonbondedForce> or <CustomGBForce> tag:

<Function name="myfn" type="Continuous1D" min="-5" max="5">
0.983674857694 -0.980096396266 -0.975743130031 -0.970451936613 -0.964027580076
-0.956237458128 -0.946806012846 -0.935409070603 -0.921668554406 -0.905148253645
-0.885351648202 -0.861723159313 -0.833654607012 -0.800499021761 -0.761594155956
-0.716297870199 -0.664036770268 -0.604367777117 -0.537049566998 -0.46211715726
-0.379948962255 -0.291312612452 -0.197375320225 -0.099667994625 0.0
0.099667994625 0.197375320225 0.291312612452 0.379948962255 0.46211715726
0.537049566998 0.604367777117 0.664036770268 0.716297870199 0.761594155956
0.800499021761 0.833654607012 0.861723159313 0.885351648202 0.905148253645
0.921668554406 0.935409070603 0.946806012846 0.956237458128 0.964027580076
0.970451936613 0.975743130031 0.980096396266 0.983674857694 0.986614298151
0.989027402201
</Function>

The tag’s attributes define the name of the function, the type of function, and the range of values for which it is defined. The required set of attributed depends on the function type:

Type

Required Attributes

Continuous1D

min, max

Continuous2D

xmin, ymin, xmax, ymax, xsize, ysize

Continuous3D

xmin, ymin, zmin, xmax, ymax, zmax, xsize, ysize, zsize

Discrete1D

Discrete2D

xsize, ysize

Discrete3D

xsize, ysize, zsize

The “min” and “max” attributes define the range of the independent variables for a continuous function. The “size” attributes define the size of the table along each axis. The tabulated values are listed inside the body of the tag, with successive values separated by white space. See the API documentation for more details.

6.2.21. Residue Template Parameters

In forces that use an <Atom> tag to define parameters for atom types or classes, there is an alternate mechanism you can also use: defining those parameter values in the residue template. This is useful for situations that come up in certain force fields. For example, NonbondedForce and GBSAOBCForce each have a charge attribute. If you only have to define the charge of each atom type once, that is more convenient and avoids potential bugs. Also, many force fields have a different charge for each atom type, but Lennard-Jones parameters that are the same for all types in a class. It would be preferable not to have to repeat those parameter values many times over.

When writing a residue template, you can add arbitrary additional attributes to each <Atom> tag. For example, you might include a charge attribute as follows:

<Atom name="CA" type="53" charge="0.0381"/>

When writing the tag for a force, you can then include a <UseAttributeFromResidue> tag inside it. This indicates that a specified attribute should be taken from the residue template. Finally, you simply omit that attribute in the force’s own <Atom> tags. For example:

<NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
 <UseAttributeFromResidue name="charge"/>
 <Atom class="CX" sigma="0.339966950842" epsilon="0.4577296"/>
 ...
</NonbondedForce>

Notice that the charge attribute is missing, and that the parameters are specified by class, not by type. This means that sigma and epsilon only need to be specified once for each class. The atom charges, which are different for each type, are taken from the residue template instead.

6.2.22. Including Other Files

Sometimes it is useful to split a force field definition into multiple files, but still be able to use the force field by specifying only a single file. You can accomplish this with the <Include> tag. For example:

<ForceField>
 <Include file="definitions.xml"/>
 ...
</ForceField>

The file attribute gives the path of the file to include. It may be relative either to the directory containing the parent XML file (the one with the <Include> tag) or the OpenMM data directory (the one containing built in force fields).

6.3. Using Multiple Files

If multiple XML files are specified when a ForceField is created, their definitions are combined as follows.

  • A file may refer to atom types and classes that it defines, as well as those defined in previous files. It may not refer to ones defined in later files. This means that the order in which files are listed when calling the ForceField constructor is potentially significant.

  • Forces that involve per-atom parameters (such as NonbondedForce or GBSAOBCForce) require parameter values to be defined for every atom type. It does not matter which file those types are defined in. For example, files that define explicit water models generally define a small number of atom types, as well as nonbonded parameters for those types. In contrast, files that define implicit solvent models do not define any new atom types, but provide parameters for all the atom types that were defined in the main force field file.

  • For other forces, the files are effectively independent. For example, if two files each include a <HarmonicBondForce> tag, bonds will be created based on the rules in the first file, and then more bonds will be created based on the rules in the second file. This means you could potentially end up with multiple bonds between a single pair of atoms.

6.4. Extending ForceField

The ForceField class is designed to be modular and extensible. This means you can add support for entirely new force types, such as ones implemented with plugins.

6.4.1. Adding new force types

For every force class, there is a “generator” class that parses the corresponding XML tag, then creates Force objects and adds them to the System. ForceField maintains a map of tag names to generator classes. When a ForceField is created, it scans through the XML files, looks up the generator class for each tag, and asks that class to create a generator object based on it. Then, when you call createSystem(), it loops over each of its generators and asks each one to create its Force object. Adding a new Force type therefore is simply a matter of creating a new generator class and adding it to ForceField’s map.

The generator class must define two methods. First, it needs a static method with the following signature to parse the XML tag and create the generator:

@staticmethod
def parseElement(element, forcefield):

element is the XML tag (an xml.etree.ElementTree.Element object) and forcefield is the ForceField being created. This method should create a generator and add it to the ForceField:

generator = MyForceGenerator()
forcefield._forces.append(generator)

It then should parse the information contained in the XML tag and configure the generator based on it.

Second, it must define a method with the following signature:

def createForce(self, system, data, nonbondedMethod, nonbondedCutoff, args):

When createSystem() is called on the ForceField, it first creates the System object, then loops over each of its generators and calls createForce() on each one. This method should create the Force object and add it to the System. data is a ForceField._SystemData object containing information about the System being created (atom types, bonds, angles, etc.), system is the System object, and the remaining arguments are values that were passed to createSystem(). To get a better idea of how this works, look at the existing generator classes in forcefield.py.

The generator class may optionally also define a method with the following signature:

def postprocessSystem(self, system, data, args):

If this method exists, it will be called after all Forces have been created. This gives generators a chance to make additional changes to the System.

Finally, you need to register your class by adding it to ForceField’s map:

forcefield.parsers['MyForce'] = MyForceGenerator.parseElement

The key is the XML tag name, and the value is the static method to use for parsing it.

Now you can simply create a ForceField object as usual. If an XML file contains a <MyForce> tag, it will be recognized and processed correctly.

6.4.2. Adding residue template generators

Caution

This feature is experimental, and its API is subject to change.

Typically, when ForceField encounters a residue it does not have a template for, it simply raises an Exception, since it does not know how to assign atom types for the unknown residue.

However, ForceField has an API for registering residue template generators that are called when a residue without an existing template is encountered. These generators may create new residue templates that match existing atom types and parameters, or can even create new atom types and new parameters that are added to ForceField. This functionality can be useful for adding residue template generators that are able to parameterize small molecules that are not represented in a protein or nucleic acid forcefield, for example, or for creating new residue templates for post-translationally modified residues, covalently-bound ligands, or unnatural amino acids or bases.

To register a new residue template generator named generator, simply call the registerTemplateGenerator() method on an existing ForceField object:

forcefield.registerTemplateGenerator(generator)

This generator function must conform to the following API:

def generator(forcefield, residue):
    """
    Parameters
    ----------
    forcefield : openmm.app.ForceField
        The ForceField object to which residue templates and/or parameters are to be added.
    residue : openmm.app.Topology.Residue
        The residue topology for which a template is to be generated.

    Returns
    -------
    success : bool
        If the generator is able to successfully parameterize the residue, `True` is returned.
        If the generator cannot parameterize the residue, it should return `False` and not
        modify `forcefield`.

    The generator should either register a residue template directly with
    `forcefield.registerResidueTemplate(template)` or it should call `forcefield.loadFile(file)`
    to load residue definitions from an ffxml file.

    It can also use the `ForceField` programmatic API to add additional atom types (via
    `forcefield.registerAtomType(parameters)`) or additional parameters.

    """

The ForceField object will be modified by the residue template generator as residues without previously defined templates are encountered. Because these templates are added to the ForceField as new residue types are encountered, subsequent residues will be parameterized using the same residue templates without calling the generator again.