GromacsTopFile supports three particle average vsites (#3634)

wrappers/python/openmm/app/gromacstopfile.py

@@ -118,6 +118,7 @@ class GromacsTopFile(object):
             self.constraints = []
             self.cmaps = []
             self.vsites2 = []
+            self.vsites3 = []
             self.has_virtual_sites = False
             self.has_nbfix_terms = False
 
@@ -257,9 +258,11 @@ class GromacsTopFile(object):
                 self._processCmapType(line)
             elif self._currentCategory == 'nonbond_params':
                 self._processNonbondType(line)
-            elif self._currentCategory == 'virtual_sites2':
+            elif self._currentCategory == 'virtual_sites2' or self._currentCategory == 'dummies2':
                 self._processVirtualSites2(line)
-            elif self._currentCategory.startswith('virtual_sites'):
+            elif self._currentCategory == 'virtual_sites3' or self._currentCategory == 'dummies3':
+                self._processVirtualSites3(line)
+            elif self._currentCategory.startswith('virtual_sites') or self._currentCategory.startswith('dummies'):
                 if self._currentMoleculeType is None:
                     raise ValueError('Found %s before [ moleculetype ]' %
                                      self._currentCategory)
@@ -468,8 +471,19 @@ class GromacsTopFile(object):
         fields = line.split()
         if len(fields) < 5:
             raise ValueError('Too few fields in [ virtual_sites2 ] line: ' + line)
+        if fields[3] != '1':
+            raise ValueError('Unsupported function type in [ virtual_sites2 ] line: '+line)
         self._currentMoleculeType.vsites2.append(fields[:5])
 
+    def _processVirtualSites3(self, line):
+        """Process a line in the [ virtual_sites3 ] category."""
+        fields = line.split()
+        if len(fields) < 7:
+            raise ValueError('Too few fields in [ virtual_sites3 ] line: ' + line)
+        if fields[4] != '1':
+            raise ValueError('Unsupported function type in [ virtual_sites3 ] line: '+line)
+        self._currentMoleculeType.vsites3.append(fields)
+
     def __init__(self, file, periodicBoxVectors=None, unitCellDimensions=None, includeDir=None, defines=None):
         """Load a top file.
 
@@ -1007,6 +1021,12 @@ class GromacsTopFile(object):
                     c1 = float(fields[4])
                     vsite = mm.TwoParticleAverageSite(baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], (1-c1), c1)
                     sys.setVirtualSite(baseAtomIndex+atoms[0], vsite)
+                for fields in moleculeType.vsites3:
+                    atoms = [int(x)-1 for x in fields[:4]]
+                    c1 = float(fields[5])
+                    c2 = float(fields[6])
+                    vsite = mm.ThreeParticleAverageSite(baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], 1-c1-c2, c1, c2)
+                    sys.setVirtualSite(baseAtomIndex+atoms[0], vsite)
 
                 # Add explicitly specified constraints.
 

wrappers/python/tests/TestGromacsTopFile.py

@@ -213,6 +213,20 @@ class TestGromacsTopFile(unittest.TestCase):
         # the energy output is from gromacs and it only prints out 6 sig digits.
         self.assertAlmostEqual(ene.value_in_unit(kilojoules_per_mole), 1.88855e+02, places=3)
 
+    def test_Vsite3(self):
+        """Test a three particle virtual site."""
+        top = GromacsTopFile('systems/tip4pew.top')
+        system = top.createSystem()
+        self.assertTrue(system.isVirtualSite(3))
+        vs = system.getVirtualSite(3)
+        self.assertIsInstance(vs, ThreeParticleAverageSite)
+        self.assertEqual(0, vs.getParticle(0))
+        self.assertEqual(1, vs.getParticle(1))
+        self.assertEqual(2, vs.getParticle(2))
+        self.assertAlmostEqual(0.786646558, vs.getWeight(0))
+        self.assertAlmostEqual(0.106676721, vs.getWeight(1))
+        self.assertAlmostEqual(0.106676721, vs.getWeight(2))
+
 if __name__ == '__main__':
     unittest.main()
 

wrappers/python/tests/systems/tip4pew.top

+; Self-contained top file representing a single TIP4P-Ew water molecule
+
+[ defaults ]
+; nbfunc	comb-rule	gen-pairs	fudgeLJ	fudgeQQ
+1		3		yes		0.5	0.5
+
+[ atomtypes ]
+; name  bond_type    mass    charge   ptype          sigma      epsilon
+MW         MW      0	0.00000     0.000       D    0.00000e+00  0.00000e+00
+HW_tip4pew   1       1.008   0.0000  A   0.00000e+00  0.00000e+00
+OW_tip4pew   8      16.00    0.0000  A   3.16435e-01  6.80946e-01
+
+[ moleculetype ]
+; molname	nrexcl
+SOL		2
+
+[ atoms ]
+; id  at type     res nr  res name  at name  cg nr  charge    mass
+  1   OW_tip4pew  1       SOL       OW       1       0        16.00000
+  2   HW_tip4pew  1       SOL       HW1      1       0.52422   1.00800
+  3   HW_tip4pew  1       SOL       HW2      1       0.52422   1.00800
+  4   MW          1       SOL       MW       1      -1.04844   0.00000
+
+[ settles ]
+; i	funct	doh	dhh
+1	1	0.09572	0.15139
+
+[ virtual_sites3 ]
+; Vsite from                    funct   a               b
+4       1       2       3       1       0.106676721     0.106676721
+
+[ exclusions ]
+1	2	3	4
+2	1	3	4
+3	1	2	4
+4	1	2	3
+
+[ molecules ]
+; Compound        #mols
+SOL     1