Add a test for 'smirnoff' impropers mode

wrappers/python/tests/TestForceField.py

@@ -907,6 +907,72 @@ 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>
+  <NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
+    <UseAttributeFromResidue name="charge"/>
+    <Atom sigma="0.3399669508423535" epsilon="0.359824" class="[H]C(=O)[H]$C1#0"/>
+    <Atom sigma="0.2959921901149463" epsilon="0.87864" class="[H]C(=O)[H]$O1#1"/>
+    <Atom sigma="0.2510552587719476" epsilon="0.06276" class="[H]C(=O)[H]$H1#2"/>
+    <Atom sigma="0.2510552587719476" epsilon="0.06276" class="[H]C(=O)[H]$H2#3"/>
+  </NonbondedForce>
+  <HarmonicBondForce>
+    <Bond class1="[H]C(=O)[H]$C1#0" class2="[H]C(=O)[H]$O1#1" length="0.122768286085" k="475133.08130977117"/>
+    <Bond class1="[H]C(=O)[H]$C1#0" class2="[H]C(=O)[H]$H1#2" length="0.10850661293899999" k="338125.5447433327"/>
+    <Bond class1="[H]C(=O)[H]$C1#0" class2="[H]C(=O)[H]$H2#3" length="0.10850661293899999" k="338125.5447433327"/>
+  </HarmonicBondForce>
+  <HarmonicAngleForce>
+    <Angle class1="[H]C(=O)[H]$O1#1" class2="[H]C(=O)[H]$C1#0" class3="[H]C(=O)[H]$H1#2" angle="2.4141195890059883" k="279.7597371776727"/>
+    <Angle class1="[H]C(=O)[H]$O1#1" class2="[H]C(=O)[H]$C1#0" class3="[H]C(=O)[H]$H2#3" angle="2.4141195890059883" k="279.7597371776727"/>
+    <Angle class1="[H]C(=O)[H]$H1#2" class2="[H]C(=O)[H]$C1#0" class3="[H]C(=O)[H]$H2#3" angle="2.2670423843424943" k="330.7820750751862"/>
+  </HarmonicAngleForce>
+  <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')