Remove support for Gromacs implicit solvent (#4246)

* Remove support for Gromacs implicit solvent

* Removed tests of implicit solvent

devtools/ci/gh-actions/conda-envs/build-macos-latest.yml

@@ -16,4 +16,4 @@ dependencies:
 - pytest
 - pytest-xdist
 - pytest-timeout
-- gromacs 2018.*
+- gromacs

devtools/ci/gh-actions/conda-envs/build-ubuntu-latest.yml

@@ -19,4 +19,4 @@ dependencies:
 - pytest
 - pytest-xdist
 - pytest-timeout
-- gromacs 2018.*
+- gromacs

devtools/ci/gh-actions/conda-envs/build-ubuntu-pypy.yml

@@ -19,4 +19,4 @@ dependencies:
 - pytest
 - pytest-xdist
 - pytest-timeout
-- gromacs 2018.*
+- gromacs

wrappers/python/openmm/app/gromacstopfile.py

@@ -290,8 +290,6 @@ class GromacsTopFile(object):
                 self._processAngleType(line)
             elif self._currentCategory == 'dihedraltypes':
                 self._processDihedralType(line)
-            elif self._currentCategory == 'implicit_genborn_params':
-                self._processImplicitType(line)
             elif self._currentCategory == 'pairtypes':
                 self._processPairType(line)
             elif self._currentCategory == 'cmaptypes':
@@ -488,13 +486,6 @@ class GromacsTopFile(object):
         else:
             self._dihedralTypes[key] = [fields]
 
-    def _processImplicitType(self, line):
-        """Process a line in the [ implicit_genborn_params ] category."""
-        fields = line.split()
-        if len(fields) < 6:
-            raise ValueError('Too few fields in [ implicit_genborn_params ] line: '+line)
-        self._implicitTypes[fields[0]] = fields
-
     def _processPairType(self, line):
         """Process a line in the [ pairtypes ] category."""
         fields = line.split()
@@ -586,7 +577,6 @@ class GromacsTopFile(object):
         self._bondTypes= {}
         self._angleTypes = {}
         self._dihedralTypes = {}
-        self._implicitTypes = {}
         self._pairTypes = {}
         self._cmapTypes = {}
         self._nonbondTypes = {}
@@ -661,9 +651,8 @@ class GromacsTopFile(object):
                 for fields in moleculeType.bonds:
                     top.addBond(atoms[int(fields[0])-1], atoms[int(fields[1])-1])
 
-    def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
-                     constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5,
-                     ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None, switchDistance=None):
+    def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer, constraints=None,
+                     rigidWater=True, ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None, switchDistance=None):
         """Construct an OpenMM System representing the topology described by this
         top file.
 
@@ -682,16 +671,6 @@ class GromacsTopFile(object):
         rigidWater : boolean=True
             If true, water molecules will be fully rigid regardless of the value
             passed for the constraints argument
-        implicitSolvent : object=None
-            If not None, the implicit solvent model to use.  The only allowed
-            value is OBC2.  This option is deprecated, since Gromacs 2019 and later
-            no longer support implicit solvent.  It will be removed in a future
-            release.
-        soluteDielectric : float=1.0
-            The solute dielectric constant to use in the implicit solvent model.
-        solventDielectric : float=78.5
-            The solvent dielectric constant to use in the implicit solvent
-            model.
         ewaldErrorTolerance : float=0.0005
             The error tolerance to use if nonbondedMethod is Ewald, PME or LJPME.
         removeCMMotion : boolean=True
@@ -741,14 +720,6 @@ class GromacsTopFile(object):
             lj.addPerParticleParameter('C')
             lj.addPerParticleParameter('A')
             sys.addForce(lj)
-        if implicitSolvent is OBC2:
-            gb = mm.GBSAOBCForce()
-            gb.setSoluteDielectric(soluteDielectric)
-            gb.setSolventDielectric(solventDielectric)
-            sys.addForce(gb)
-            nb.setReactionFieldDielectric(1.0)
-        elif implicitSolvent is not None:
-            raise ValueError('Illegal value for implicitSolvent')
         bonds = {}
         angles = {}
         periodic = None
@@ -1038,11 +1009,6 @@ class GromacsTopFile(object):
                             epsilon = float(params[7])
                             lj.addParticle([math.sqrt(4*epsilon*sigma**6), math.sqrt(4*epsilon*sigma**12)])
 
-                    if implicitSolvent is OBC2:
-                        if fields[1] not in self._implicitTypes:
-                            raise ValueError('No implicit solvent parameters specified for atom type: '+fields[1])
-                        gbparams = self._implicitTypes[fields[1]]
-                        gb.addParticle(q, float(gbparams[4]), float(gbparams[5]))
                 for fields in moleculeType.bonds:
                     atoms = [int(x)-1 for x in fields[:2]]
                     bondIndices.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1]))

wrappers/python/tests/TestGromacsTopFile.py

@@ -19,8 +19,6 @@ class TestGromacsTopFile(unittest.TestCase):
 
         # alanine dipeptide with explicit water
         self.top1 = GromacsTopFile('systems/explicit.top', unitCellDimensions=Vec3(6.223, 6.223, 6.223)*nanometers)
-        # alanine dipeptide with implicit water
-        self.top2 = GromacsTopFile('systems/implicit.top')
 
     def test_NonbondedMethod(self):
         """Test all six options for the nonbondedMethod parameter."""
@@ -96,17 +94,16 @@ class TestGromacsTopFile(unittest.TestCase):
 
     def test_SwitchingFunction(self):
         """Test using a switching function."""
-        for filename in ('systems/implicit.top', 'systems/ionic.top'):
-            top = GromacsTopFile(filename)
-            for distance in (None, 0.8*nanometers):
-                system = top.createSystem(nonbondedMethod=CutoffNonPeriodic, switchDistance=distance)
-                for f in system.getForces():
-                    if isinstance(f, NonbondedForce) or isinstance(f, CustomNonbondedForce):
-                        if distance is None:
-                            self.assertFalse(f.getUseSwitchingFunction())
-                        else:
-                            self.assertTrue(f.getUseSwitchingFunction())
-                            self.assertEqual(distance, f.getSwitchingDistance())
+        top = GromacsTopFile('systems/ionic.top')
+        for distance in (None, 0.8*nanometers):
+            system = top.createSystem(nonbondedMethod=CutoffNonPeriodic, switchDistance=distance)
+            for f in system.getForces():
+                if isinstance(f, NonbondedForce) or isinstance(f, CustomNonbondedForce):
+                    if distance is None:
+                        self.assertFalse(f.getUseSwitchingFunction())
+                    else:
+                        self.assertTrue(f.getUseSwitchingFunction())
+                        self.assertEqual(distance, f.getSwitchingDistance())
 
     def test_EwaldErrorTolerance(self):
         """Test to make sure the ewaldErrorTolerance parameter is passed correctly."""
@@ -140,37 +137,6 @@ class TestGromacsTopFile(unittest.TestCase):
                 validateConstraints(self, topology, system,
                                     constraints_value, rigidWater_value)
 
-    def test_ImplicitSolvent(self):
-        """Test implicit solvent using the implicitSolvent parameter.
-
-        """
-        system = self.top2.createSystem(implicitSolvent=OBC2)
-        self.assertTrue(any(isinstance(f, GBSAOBCForce) for f in system.getForces()))
-
-    def test_ImplicitSolventParameters(self):
-        """Test that solventDielectric and soluteDielectric are passed correctly.
-
-        """
-        system = self.top2.createSystem(implicitSolvent=OBC2,
-                                           solventDielectric=50.0,
-                                           soluteDielectric = 0.9)
-        found_matching_solvent_dielectric=False
-        found_matching_solute_dielectric=False
-        for force in system.getForces():
-            if isinstance(force, GBSAOBCForce):
-                if force.getSolventDielectric() == 50.0:
-                    found_matching_solvent_dielectric = True
-                if force.getSoluteDielectric() == 0.9:
-                    found_matching_solute_dielectric = True
-                gbcharges = [force.getParticleParameters(i)[0] for i in range(system.getNumParticles())]
-            if isinstance(force, NonbondedForce):
-                self.assertEqual(force.getReactionFieldDielectric(), 1.0)
-                nbcharges = [force.getParticleParameters(i)[0] for i in range(system.getNumParticles())]
-        self.assertTrue(found_matching_solvent_dielectric and
-                        found_matching_solute_dielectric)
-        for q1, q2 in zip(gbcharges, nbcharges):
-            self.assertEqual(q1, q2)
-
     def test_HydrogenMass(self):
         """Test that altering the mass of hydrogens works correctly."""