Merge pull request #2511 from jchodera/smirnoff-impropers

Add support to ForceField for handling SMIRNOFF trefoil impropers

docs-source/usersguide/application.rst

@@ -2331,6 +2331,22 @@ 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 :code:`<PeriodicTorsionForce>` tag also supports an optional
+:code:`ordering` attribute to provide better compatibility with the way
+impropers are assigned in different simulation packages:
+
+ * :code:`ordering="default"` specifies the default behavior if the attribute
+   is omitted. 
+ * :code:`ordering="amber"` produces behavior that replicates the behavior of
+   AmberTools LEaP
+ * :code:`ordering="charmm"` produces behavior more consistent with CHARMM
+ * :code:`ordering="smirnoff"` allows multiple impropers to be added using
+   exact matching to replicate the beheavior of `SMIRNOFF <https://open-forcefield-toolkit.readthedocs.io/en/latest/smirnoff.html>`_
+   impropers
+
+Different :code:`<PeriodicTorsionForce>` tags can specify different :code:`ordering`
+values to be used for the sub-elements appearing within their tags.
+
 <RBTorsionForce>
 ================
 

wrappers/python/simtk/openmm/app/forcefield.py

@@ -1841,6 +1841,16 @@ def _matchImproper(data, torsion, generator):
                                 (a2, a3) = (a3, a2)
                         match = (a2, a3, torsion[0], a4, tordef)
                         break
+                    elif tordef.ordering == 'smirnoff':
+                        # topology atom indexes
+                        a1 = torsion[0]
+                        a2 = torsion[t2[1]]
+                        a3 = torsion[t3[1]]
+                        a4 = torsion[t4[1]]
+                        # enforce exact match
+                        match = (a1, a2, a3, a4, tordef)
+                        break
+
     return match
 
 
@@ -2032,7 +2042,7 @@ class PeriodicTorsionGenerator(object):
     def registerImproperTorsion(self, parameters, ordering='default'):
         torsion = self.ff._parseTorsion(parameters)
         if torsion is not None:
-            if ordering in ['default', 'charmm', 'amber']:
+            if ordering in ['default', 'charmm', 'amber', 'smirnoff']:
                 torsion.ordering = ordering
             else:
                 raise ValueError('Illegal ordering type %s for improper torsion %s' % (ordering, torsion))
@@ -2116,6 +2126,12 @@ class PeriodicTorsionGenerator(object):
                 (a1, a2, a3, a4, tordef) = match
                 for i in range(len(tordef.phase)):
                     if tordef.k[i] != 0:
+                        if tordef.ordering == 'smirnoff':
+                            # Add all torsions in trefoil
+                            force.addTorsion(a1, a2, a3, a4, tordef.periodicity[i], tordef.phase[i], tordef.k[i])
+                            force.addTorsion(a1, a3, a4, a2, tordef.periodicity[i], tordef.phase[i], tordef.k[i])
+                            force.addTorsion(a1, a4, a2, a3, tordef.periodicity[i], tordef.phase[i], tordef.k[i])
+                        else:
                             force.addTorsion(a1, a2, a3, a4, tordef.periodicity[i], tordef.phase[i], tordef.k[i])
 parsers["PeriodicTorsionForce"] = PeriodicTorsionGenerator.parseElement
 

wrappers/python/tests/TestForceField.py

@@ -907,6 +907,55 @@ class TestForceField(unittest.TestCase):
         self.assertEqual(system1_indexes, [51, 56, 54, 55])
         self.assertEqual(system2_indexes, [51, 55, 54, 56])
 
+    def test_ImpropersOrdering_smirnoff(self):
+        """Test correctness of the ordering of atom indexes in improper torsions
+        and the torsion.ordering parameter when using the 'smirnoff' mode.
+        """
+
+        # SMIRNOFF parameters for formaldehyde
+        xml = """
+<ForceField>
+  <AtomTypes>
+    <Type name="[H]C(=O)[H]$C1#0" element="C" mass="12.01078" class="[H]C(=O)[H]$C1#0"/>
+    <Type name="[H]C(=O)[H]$O1#1" element="O" mass="15.99943" class="[H]C(=O)[H]$O1#1"/>
+    <Type name="[H]C(=O)[H]$H1#2" element="H" mass="1.007947" class="[H]C(=O)[H]$H1#2"/>
+    <Type name="[H]C(=O)[H]$H2#3" element="H" mass="1.007947" class="[H]C(=O)[H]$H2#3"/>
+  </AtomTypes>
+  <PeriodicTorsionForce ordering="smirnoff">
+    <Improper class1="[H]C(=O)[H]$C1#0" class2="[H]C(=O)[H]$O1#1" class3="[H]C(=O)[H]$H1#2" class4="[H]C(=O)[H]$H2#3" periodicity1="2" phase1="3.141592653589793" k1="1.5341333333333336"/>
+    <Improper class1="[H]C(=O)[H]$C1#0" class2="[H]C(=O)[H]$H1#2" class3="[H]C(=O)[H]$H2#3" class4="[H]C(=O)[H]$O1#1" periodicity1="2" phase1="3.141592653589793" k1="1.5341333333333336"/>
+    <Improper class1="[H]C(=O)[H]$C1#0" class2="[H]C(=O)[H]$H2#3" class3="[H]C(=O)[H]$O1#1" class4="[H]C(=O)[H]$H1#2" periodicity1="2" phase1="3.141592653589793" k1="1.5341333333333336"/>
+  </PeriodicTorsionForce>
+  <Residues>
+    <Residue name="[H]C(=O)[H]">
+      <Atom name="C1" type="[H]C(=O)[H]$C1#0" charge="0.5632799863815308"/>
+      <Atom name="O1" type="[H]C(=O)[H]$O1#1" charge="-0.514739990234375"/>
+      <Atom name="H1" type="[H]C(=O)[H]$H1#2" charge="-0.02426999807357788"/>
+      <Atom name="H2" type="[H]C(=O)[H]$H2#3" charge="-0.02426999807357788"/>
+      <Bond atomName1="C1" atomName2="O1"/>
+      <Bond atomName1="C1" atomName2="H1"/>
+      <Bond atomName1="C1" atomName2="H2"/>
+    </Residue>
+  </Residues>
+</ForceField>
+"""
+        pdb = PDBFile('systems/formaldehyde.pdb')
+        # ff1 uses default ordering of impropers, ff2 uses "amber" for the one
+        # problematic improper
+        ff = ForceField(StringIO(xml))
+
+        system = ff.createSystem(pdb.topology)
+
+        # Check that impropers are applied in the correct three-fold trefoil pattern
+        forces = { force.__class__.__name__ : force for force in system.getForces() }
+        force = forces['PeriodicTorsionForce']
+        created_torsions = set()
+        for index in range(force.getNumTorsions()):
+            i,j,k,l,_,_,_ = force.getTorsionParameters(index)
+            created_torsions.add((i,j,k,l))
+        expected_torsions = set([(0,3,1,2), (0,1,2,3), (0,2,3,1)])
+        self.assertEqual(expected_torsions, created_torsions)
+
     def test_Disulfides(self):
         """Test that various force fields handle disulfides correctly."""
         pdb = PDBFile('systems/bpti.pdb')

wrappers/python/tests/systems/formaldehyde.pdb

+REMARK   1 CREATED WITH OPENMM 7.4, 2020-01-04
+HETATM    1          A   1       0.095   0.011   0.000  1.00  0.00           C  
+HETATM    2          A   1       0.513  -1.098  -0.004  1.00  0.00           O  
+HETATM    3          A   1      -1.098   0.149   0.015  1.00  0.00           H  
+HETATM    4          A   1       0.590   0.938  -0.011  1.00  0.00           H  
+TER       5          A   1
+CONECT    1    2    3    4
+CONECT    2    1
+CONECT    3    1
+CONECT    4    1
+END