GromacsTopFile supports vsite type 3fad (#5155)

wrappers/python/openmm/app/gromacstopfile.py

@@ -4,7 +4,7 @@ gromacstopfile.py: Used for loading Gromacs top files.
 This is part of the OpenMM molecular simulation toolkit.
 See https://openmm.org/development.
 
-Portions copyright (c) 2012-2024 Stanford University and the Authors.
+Portions copyright (c) 2012-2025 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors: Jason Swails
 
@@ -525,7 +525,7 @@ class GromacsTopFile(object):
         fields = line.split()
         if len(fields) < 7:
             raise ValueError('Too few fields in [ virtual_sites3 ] line: ' + line)
-        if fields[4] not in ('1', '4'):
+        if fields[4] not in ('1', '3', '4'):
             raise ValueError('Unsupported function type in [ virtual_sites3 ] line: '+line)
         self._currentMoleculeType.vsites3.append(fields)
 
@@ -1117,6 +1117,10 @@ class GromacsTopFile(object):
                     c2 = float(fields[6])
                     if vsiteType == '1':
                         vsite = mm.ThreeParticleAverageSite(baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], 1-c1-c2, c1, c2)
+                    elif vsiteType == '3':
+                        theta = c1*math.pi/180
+                        vsite = mm.LocalCoordinatesSite(baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3],
+                                                        (1, 0, 0), (-1, 1, 0), (0, -1, 1), (c2*math.cos(theta), c2*math.sin(theta), 0))
                     elif vsiteType == '4':
                         c3 = float(fields[7])
                         vsite = mm.OutOfPlaneSite(baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], c1, c2, c3)

wrappers/python/tests/TestGromacsTopFile.py

@@ -225,6 +225,18 @@ class TestGromacsTopFile(unittest.TestCase):
         self.assertAlmostEqual(0.106676721, vs.getWeight(1))
         self.assertAlmostEqual(0.106676721, vs.getWeight(2))
 
+    def test_Vsite3Func3(self):
+        """Test a three particle virtual site with function 3."""
+        top = GromacsTopFile('systems/vsite3-3.top')
+        gro = GromacsGroFile('systems/vsite3-3.gro')
+        system = top.createSystem()
+        context = Context(system, VerletIntegrator(1.0), Platform.getPlatform('Reference'))
+        context.setPositions(gro.positions)
+        context.computeVirtualSites()
+        positions = context.getState(positions=True).getPositions().value_in_unit(nanometer)
+        # Compare to the position computed by Gromacs
+        assert_allclose([2.45, 0.7794, 0.0], positions[3], atol=1e-4)
+
     def test_Vsite3Func4(self):
         """Test a three particle virtual site."""
         top = GromacsTopFile('systems/tip5p.top')

wrappers/python/tests/systems/vsite3-3.gro

+Vsite test system
+    4
+    1MOL     C1    1   0.000   1.000   0.000
+    1MOL     C2    2   1.000   0.000   0.000
+    1MOL     C3    3   2.000   0.000   0.000
+    1MOL     H     4   3.000   0.000   0.000
+   8.00000   8.00000   8.00000

wrappers/python/tests/systems/vsite3-3.top

+[ defaults ]
+; nbfunc	comb-rule	gen-pairs	fudgeLJ	fudgeQQ
+1		3		yes		0.5	0.5
+
+[ atomtypes ]
+; name  bond_type    mass    charge   ptype          sigma      epsilon
+H   1       0.00    0.0000  A   0.00000e+00  0.00000e+00
+C   6      12.00    0.0000  A   0.00000e+00  0.00000e+00
+
+[ moleculetype ]
+; molname	nrexcl
+MOL		3
+
+[ atoms ]
+; id  at type     res nr  res name  at name  cg nr  charge    mass
+  1   C           1       MOL       C1       1       0        12.00000
+  2   C           1       MOL       C2       1       0        12.00000
+  3   C           1       MOL       C3       1       0        12.00000
+  4   H           1       MOL       H        1       0         0.00000
+
+[ bonds ]
+
+[ virtual_sites3 ]
+; Vsite from                    funct   theta   d
+4       3       2       1       3       120     0.9
+
+[ exclusions ]
+1	2	3	4
+2	1	3	4
+3	1	2	4
+4	1	2	3
+
+[ system ]
+1 MOL
+
+[ molecules ]
+; Compound        #mols
+MOL     1