GromacsTopFile supports virtual_sites3 function 4 (#4536)

wrappers/python/openmm/app/gromacstopfile.py

@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
 Biological Structures at Stanford, funded under the NIH Roadmap for
 Medical Research, grant U54 GM072970. See https://simtk.org.
 
-Portions copyright (c) 2012-2023 Stanford University and the Authors.
+Portions copyright (c) 2012-2024 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors: Jason Swails
 
@@ -527,7 +527,7 @@ class GromacsTopFile(object):
         fields = line.split()
         if len(fields) < 7:
             raise ValueError('Too few fields in [ virtual_sites3 ] line: ' + line)
-        if fields[4] != '1':
+        if fields[4] not in ('1', '4'):
             raise ValueError('Unsupported function type in [ virtual_sites3 ] line: '+line)
         self._currentMoleculeType.vsites3.append(fields)
 
@@ -1079,9 +1079,16 @@ class GromacsTopFile(object):
                     sys.setVirtualSite(baseAtomIndex+atoms[0], vsite)
                 for fields in moleculeType.vsites3:
                     atoms = [int(x)-1 for x in fields[:4]]
+                    vsiteType = 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)
+                    if vsiteType == '1':
+                        vsite = mm.ThreeParticleAverageSite(baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], 1-c1-c2, c1, c2)
+                    elif vsiteType == '4':
+                        c3 = float(fields[7])
+                        vsite = mm.OutOfPlaneSite(baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], c1, c2, c3)
+                    else:
+                        raise ValueError('Internal error: vsites3 has unexpected type: '+vsiteType)
                     sys.setVirtualSite(baseAtomIndex+atoms[0], vsite)
 
                 # Add explicitly specified constraints.

wrappers/python/tests/TestGromacsTopFile.py

@@ -180,7 +180,7 @@ 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):
+    def test_Vsite3Func1(self):
         """Test a three particle virtual site."""
         top = GromacsTopFile('systems/tip4pew.top')
         system = top.createSystem()
@@ -195,6 +195,25 @@ class TestGromacsTopFile(unittest.TestCase):
         self.assertAlmostEqual(0.106676721, vs.getWeight(1))
         self.assertAlmostEqual(0.106676721, vs.getWeight(2))
 
+    def test_Vsite3Func4(self):
+        """Test a three particle virtual site."""
+        top = GromacsTopFile('systems/tip5p.top')
+        system = top.createSystem()
+        self.assertEqual(3, system.getNumConstraints())
+        for i in (3, 4):
+            self.assertTrue(system.isVirtualSite(i))
+            vs = system.getVirtualSite(i)
+            self.assertIsInstance(vs, OutOfPlaneSite)
+            self.assertEqual(0, vs.getParticle(0))
+            self.assertEqual(1, vs.getParticle(1))
+            self.assertEqual(2, vs.getParticle(2))
+            self.assertAlmostEqual(-0.344908, vs.getWeight12())
+            self.assertAlmostEqual(-0.344908, vs.getWeight13())
+            wc = -6.4437903493
+            if i == 4:
+                wc = -wc
+            self.assertAlmostEqual(wc, vs.getWeightCross())
+
     def test_GROMOS(self):
         """Test a system using the GROMOS 54a7 force field."""
 

wrappers/python/tests/systems/tip5p.top

+#include "oplsaa.ff/forcefield.itp"
+
+; original reference: [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112 , 2000]
+;
+; Note that there are various issues with tip5p and the different forcefields.
+; Discussion is here: https://gitlab.com/gromacs/gromacs/-/issues/1348
+
+[ moleculetype ]
+; molname       nrexcl
+SOL             2
+
+[ atoms ]
+; id    at type res nr  residu name     at name         cg nr   charge
+1       opls_118     1       SOL              OW             1       0
+2       opls_119     1       SOL             HW1             1       0.241
+3       opls_119     1       SOL             HW2             1       0.241
+4       opls_120     1       SOL             LP1             1      -0.241
+5       opls_120     1       SOL             LP2             1      -0.241
+
+[ bonds ]
+; i     j       funct   length  force.c.
+1       2       1       0.09572 502416.0 0.09572        502416.0
+1       3       1       0.09572 502416.0 0.09572        502416.0
+
+[ angles ]
+; i     j       k       funct   angle   force.c.
+2       1       3       1       104.52  628.02  104.52  628.02
+
+[ virtual_sites3 ]
+; The position of the virtual site is computed as follows:
+;
+; The distance from OW to L is 0.07 nm, the geometry is tetrahedral
+; (109.47 deg)
+; Therefore, a = b = 0.07 * cos (109.47/2) / | xOH1 + xOH2 |
+;            c = 0.07 * sin (109.47/2) / | xOH1 X xOH2 |
+;
+; Using | xOH1 X xOH2 | = | xOH1 | | xOH2 | sin (H1-O-H2)
+;       | xOH1 + xOH2 | = 2 | xOH1 | cos (H1-O-H2)
+; Vsite pos x4 = x1 + a*x21 + b*x31 + c*(x21 X x31)
+
+; Vsite from                    funct   a       b               c
+4       1       2       3       4       -0.344908  -0.344908  -6.4437903493
+5       1       2       3       4       -0.344908  -0.344908   6.4437903493
+
+[ exclusions ]
+1       2       3       4       5
+2       1       3       4       5
+3       1       2       4       5
+4       1       2       3       5
+5       1       2       3       4
+
+[ system ]
+1 TIP5P
+
+[ molecules ]
+SOL 1