Merge pull request #1182 from peastman/residueparams

[WIP] Per-atom parameters can be taken from residue templates

docs-source/usersguide/application.rst

@@ -1843,7 +1843,7 @@ The :code:`<ForceField>` tag contains the following children:
 * Zero or more tags defining specific forces
 
 
-The order of these tags does not matter.  They are described in details below.
+The order of these tags does not matter.  They are described in detail below.
 
 <AtomTypes>
 ===========
@@ -1927,9 +1927,9 @@ as in the following example:
 
 Each :code:`<VirtualSite>` tag indicates an atom in the residue that should
 be represented with a virtual site.  The :code:`type` attribute may equal
-:code:`"average2"`\ , :code:`"average3"`\ , or :code:`"outOfPlane"`\ , which
-correspond to the TwoParticleAverageSite, ThreeParticleAverageSite, and
-OutOfPlaneSite classes respectively.  The :code:`index` attribute gives the
+:code:`"average2"`\ , :code:`"average3"`\ , :code:`"outOfPlane"`\ , or
+:code:`"localCoords"`\ , which correspond to the TwoParticleAverageSite, ThreeParticleAverageSite,
+OutOfPlaneSite, and LocalCoordinatesSite classes respectively.  The :code:`index` attribute gives the
 index (starting from 0) of the atom to represent with a virtual site.  The atoms
 it is calculated based on are specified by :code:`atom1`\ , :code:`atom2`\ ,
 and (for virtual site classes that involve three atoms) :code:`atom3`\ .  The
@@ -1938,7 +1938,10 @@ parameters for calculating the site position.  For a TwoParticleAverageSite,
 they are :code:`weight1` and :code:`weight2`\ .  For a
 ThreeParticleAverageSite, they are :code:`weight1`\ , :code:`weight2`\ , and
 \ :code:`weight3`\ . For an OutOfPlaneSite, they are :code:`weight12`\ ,
-:code:`weight13`\ , and :code:`weightCross`\ .
+:code:`weight13`\ , and :code:`weightCross`\ . For a LocalCoordinatesSite, they
+are :code:`wo1`\ , :code:`wo2`\ , :code:`wo3`\ , :code:`wx1`\ , :code:`wx2`\ ,
+:code:`wx3`\ , :code:`wy1`\ , :code:`wy2`\ , :code:`wy3`\ , :code:`p1`\ ,
+:code:`p2`\ , and :code:`p3`\ .
 
 
 <HarmonicBondForce>
@@ -2496,8 +2499,9 @@ The following operators are supported: + (add), - (subtract), * (multiply), /
 
 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,
-step. step(x) = 0 if x < 0, 1 otherwise.  Some custom forces allow additional
-functions to be defined from tabulated values.
+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.
@@ -2521,8 +2525,8 @@ values.  All uses of a value must appear *before* that value’s definition.
 
 .. _tabulated-functions:
 
-TabulatedFunctions
-==================
+Tabulated Functions
+===================
 
 Some forces, such as CustomNonbondedForce and CustomGBForce, allow you to define
 tabulated functions.  To define a function, include a :code:`<Function>` tag inside the
@@ -2569,6 +2573,47 @@ successive values separated by white space.  See the API documentation for more
 details.
 
 
+Residue Template Parameters
+===========================
+
+In forces that use an :code:`<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, :code:`NonbondedForce` and
+:code:`GBSAOBCForce` each have a :code:`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 :code:`<Atom>` tag.  For example, you might include a :code:`charge`
+attribute as follows:
+
+.. code-block:: xml
+
+   <Atom name="CA" type="53" charge="0.0381"/>
+
+When writing the tag for a force, you can then include a
+:code:`<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 :code:`<Atom>` tags.  For example:
+
+.. code-block:: xml
+
+    <NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
+     <UseAttributeFromResidue name="charge"/>
+     <Atom class="CX" sigma="0.339966950842" epsilon="0.4577296"/>
+     ...
+    </NonbondedForce>
+
+Notice that the :code:`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.
+
+
 Using Multiple Files
 ********************
 

wrappers/python/tests/TestForceField.py

@@ -6,6 +6,10 @@ from simtk.unit import *
 import simtk.openmm.app.element as elem
 import simtk.openmm.app.forcefield as forcefield
 import math
+if sys.version_info >= (3, 0):
+    from io import StringIO
+else:
+    from cStringIO import StringIO
 
 class TestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method."""
@@ -144,7 +148,8 @@ class TestForceField(unittest.TestCase):
         context = Context(system, integrator)
         context.setPositions(pdb.positions)
         state1 = context.getState(getForces=True)
-        state2 = XmlSerializer.deserialize(open('systems/lysozyme-implicit-forces.xml').read())
+        with open('systems/lysozyme-implicit-forces.xml') as input:
+            state2 = XmlSerializer.deserialize(input.read())
         numDifferences = 0
         for f1, f2, in zip(state1.getForces().value_in_unit(kilojoules_per_mole/nanometer), state2.getForces().value_in_unit(kilojoules_per_mole/nanometer)):
             diff = norm(f1-f2)
@@ -198,6 +203,50 @@ class TestForceField(unittest.TestCase):
             self.assertEqual(vectors[i], self.pdb1.topology.getPeriodicBoxVectors()[i])
             self.assertEqual(vectors[i], system.getDefaultPeriodicBoxVectors()[i])
 
+    def test_ResidueAttributes(self):
+        """Test a ForceField that gets per-particle parameters from residue attributes."""
+
+        xml = """
+<ForceField>
+ <AtomTypes>
+  <Type name="tip3p-O" class="OW" element="O" mass="15.99943"/>
+  <Type name="tip3p-H" class="HW" element="H" mass="1.007947"/>
+ </AtomTypes>
+ <Residues>
+  <Residue name="HOH">
+   <Atom name="O" type="tip3p-O" charge="-0.834"/>
+   <Atom name="H1" type="tip3p-H" charge="0.417"/>
+   <Atom name="H2" type="tip3p-H" charge="0.417"/>
+   <Bond from="0" to="1"/>
+   <Bond from="0" to="2"/>
+  </Residue>
+ </Residues>
+ <NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
+  <UseAttributeFromResidue name="charge"/>
+  <Atom type="tip3p-O" sigma="0.315" epsilon="0.635"/>
+  <Atom type="tip3p-H" sigma="1" epsilon="0"/>
+ </NonbondedForce>
+</ForceField>"""
+        ff = ForceField(StringIO(xml))
+
+        # Build a water box.
+        modeller = Modeller(Topology(), [])
+        modeller.addSolvent(ff, boxSize=Vec3(3, 3, 3)*nanometers)
+
+        # Create a system and make sure all nonbonded parameters are correct.
+        system = ff.createSystem(modeller.topology)
+        nonbonded = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]
+        atoms = list(modeller.topology.atoms())
+        for i in range(len(atoms)):
+            params = nonbonded.getParticleParameters(i)
+            if atoms[i].element == elem.oxygen:
+                self.assertEqual(params[0], -0.834*elementary_charge)
+                self.assertEqual(params[1], 0.315*nanometers)
+                self.assertEqual(params[2], 0.635*kilojoule_per_mole)
+            else:
+                self.assertEqual(params[0], 0.417*elementary_charge)
+                self.assertEqual(params[1], 1.0*nanometers)
+                self.assertEqual(params[2], 0.0*kilojoule_per_mole)
 
 class AmoebaTestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method with the AMOEBA forcefield."""
@@ -270,6 +319,22 @@ class AmoebaTestForceField(unittest.TestCase):
         self.assertAlmostEqual(constraints[(0,2)], hoDist)
         self.assertAlmostEqual(constraints[(1,2)], hohDist)
 
+    def test_Forces(self):
+        """Compute forces and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
+
+        pdb = PDBFile('systems/alanine-dipeptide-implicit.pdb')
+        forcefield = ForceField('amoeba2013.xml', 'amoeba2013_gk.xml')
+        system = forcefield.createSystem(pdb.topology, polarization='direct')
+        integrator = VerletIntegrator(0.001)
+        context = Context(system, integrator)
+        context.setPositions(pdb.positions)
+        state1 = context.getState(getForces=True)
+        with open('systems/alanine-dipeptide-amoeba-forces.xml') as input:
+            state2 = XmlSerializer.deserialize(input.read())
+        for f1, f2, in zip(state1.getForces().value_in_unit(kilojoules_per_mole/nanometer), state2.getForces().value_in_unit(kilojoules_per_mole/nanometer)):
+            diff = norm(f1-f2)
+            self.assertTrue(diff < 0.1 or diff/norm(f1) < 1e-3)
+
 if __name__ == '__main__':
     unittest.main()
 

wrappers/python/tests/systems/alanine-dipeptide-amoeba-forces.xml

+<?xml version="1.0" ?>
+<State time="0" type="State" version="1">
+	<PeriodicBoxVectors>
+		<A x="2" y="0" z="0"/>
+		<B x="0" y="2" z="0"/>
+		<C x="0" y="0" z="2"/>
+	</PeriodicBoxVectors>
+	<Forces>
+		<Force x="-60.899006905648534" y="-704.8991281062977" z=".0364180169671772"/>
+		<Force x="1179.7263311404506" y="629.2168542831286" z="1.563821009034657"/>
+		<Force x="-484.4401097634541" y="135.94210610096957" z="558.0155391926069"/>
+		<Force x="-484.1625627429897" y="136.02915381525668" z="-557.7203977673528"/>
+		<Force x="-804.8403408476064" y="-857.579131840898" z="118.43755474210333"/>
+		<Force x="84.24441920455223" y="-421.86322192608037" z="120.734742628977"/>
+		<Force x="666.4942549112259" y="-479.03339298154526" z="19.66166555862896"/>
+		<Force x="-518.9062723002917" y="686.3850843706331" z=".7311137123231983"/>
+		<Force x="-105.62094199225953" y="710.9845370176521" z="-508.0677463920413"/>
+		<Force x="285.6690654282444" y="-111.35548608152155" z="847.8096733172159"/>
+		<Force x="-396.4822698774795" y="24.462223992557014" z="297.9581483103431"/>
+		<Force x="-168.22340166944585" y="110.99779797739818" z="-824.0916382722021"/>
+		<Force x="805.7023041833917" y="15.738288953519556" z="-340.2263402925698"/>
+		<Force x="366.3025998383222" y="-647.7535895818237" z="-273.83047817845636"/>
+		<Force x="-649.7310636705703" y="-495.52797252280067" z="317.36013856851315"/>
+		<Force x="-480.120788018345" y="172.44477988154952" z="-16.235455767965004"/>
+		<Force x="-486.20230228601736" y="982.6376366390087" z="127.20818643942287"/>
+		<Force x="820.5730252078641" y="-78.71234329038225" z="109.31395566365498"/>
+		<Force x="355.03388135028814" y="-1391.2776654673507" z=".9942543092954484"/>
+		<Force x="-511.3467906691061" y="677.7717086552518" z=".008253291392112616"/>
+		<Force x="293.66157422204157" y="452.68343054060256" z="502.83463917630775"/>
+		<Force x="293.5683952568295" y="452.70832957117113" z="-502.49604726620026"/>
+	</Forces>
+</State>