Forbid LJPME when LJ uses CustomNonbondedForce (#5044)

wrappers/python/openmm/app/charmmpsffile.py

@@ -1217,6 +1217,8 @@ class CharmmPsfFile(object):
                 force.addParticle(atm.charge, sigma_scale*atm.type.rmin*length_conv,
                                   abs(atm.type.epsilon*ene_conv))
         else:
+            if nonbondedMethod is ff.LJPME:
+                raise ValueError('LJPME is not supported when NBFIX terms are present')
             for atm in self.atom_list:
                 force.addParticle(atm.charge, 1.0, 0.0)
             # Now add the custom nonbonded force that implements NBFIX. First

wrappers/python/openmm/app/desmonddmsfile.py

@@ -368,12 +368,14 @@ class DesmondDMSFile(object):
         nb.setEwaldErrorTolerance(ewaldErrorTolerance)
         if cnb is not None:
             nb.setUseDispersionCorrection(False)
-            if nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME, ff.LJPME):
+            if nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME):
                 cnb.setNonbondedMethod(methodMap[ff.CutoffPeriodic])
                 cnb.setCutoffDistance(nonbondedCutoff)
             elif nonbondedMethod == ff.CutoffNonPeriodic:
                 cnb.setNonbondedMethod(methodMap[ff.CutoffNonPeriodic])
                 cnb.setCutoffDistance(nonbondedCutoff)
+            elif nonbondedMethod is ff.LJPME:
+                raise ValueError('LJPME is not supported with OPLS combining rules')
             else:
                 cnb.setNonbondedMethod(methodMap[ff.NoCutoff])
             cnb.setUseSwitchingFunction(False)

wrappers/python/openmm/app/forcefield.py

@@ -2836,12 +2836,14 @@ class LennardJonesGenerator(object):
         self.force.addTabulatedFunction('bcoef', mm.Discrete2DFunction(numLjTypes, numLjTypes, bcoef))
         self.force.addPerParticleParameter('type')
         self.force.setName('LennardJones')
-        if nonbondedMethod in [CutoffPeriodic, Ewald, PME, LJPME]:
+        if nonbondedMethod in [CutoffPeriodic, Ewald, PME]:
             self.force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic)
         elif nonbondedMethod is NoCutoff:
             self.force.setNonbondedMethod(mm.CustomNonbondedForce.NoCutoff)
         elif nonbondedMethod is CutoffNonPeriodic:
             self.force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic)
+        elif nonbondedMethod is LJPME:
+            raise ValueError('LJPME is not supported by LennardJonesForce')
         else:
             raise AssertionError('Unrecognized nonbonded method [%s]' % nonbondedMethod)
         if args['switchDistance'] is not None:

wrappers/python/openmm/app/gromacstopfile.py

@@ -1239,6 +1239,8 @@ class GromacsTopFile(object):
             nb.setUseSwitchingFunction(True)
             nb.setSwitchingDistance(switchDistance)
 
+        if not (lj is None and ljnbfix is None) and nonbondedMethod is ff.LJPME:
+            raise ValueError('LJPME is not supported with active nonbond_params or combining rules other than 2')
         if lj is not None:
             methodMap = {ff.NoCutoff:mm.CustomNonbondedForce.NoCutoff,
                          ff.CutoffNonPeriodic:mm.CustomNonbondedForce.CutoffNonPeriodic,

wrappers/python/openmm/app/internal/amber_file_parser.py

@@ -1037,6 +1037,8 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
 
     # Copy the exceptions as exclusions to the CustomNonbondedForce if we have
     # NBFIX terms
+    if nbfix and nonbondedMethod == 'LJPME':
+        raise ValueError('LJPME is not supported with modified off-diagonal Lennard-Jones coefficients')
     if nbfix or has_1264:
         for i in range(force.getNumExceptions()):
             ii, jj, chg, sig, eps = force.getExceptionParameters(i)

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -39,6 +39,12 @@ class TestAmberPrmtopFile(unittest.TestCase):
                                 f.getNonbondedMethod()==methodMap[method]
                                 for f in forces))
 
+    def test_LJPME_NBFIX(self):
+        """Test LJPME with NBFIX."""
+
+        with self.assertRaisesRegex(ValueError, 'LJPME is not supported'):
+            prmtop3.createSystem(nonbondedMethod=LJPME)
+
     def test_Cutoff(self):
         """Test to make sure the nonbondedCutoff parameter is passed correctly."""
 

wrappers/python/tests/TestCharmmFiles.py

@@ -287,6 +287,8 @@ class TestCharmmFiles(unittest.TestCase):
         # Check some illegal options
         self.assertRaises(ValueError, lambda:
                     psf.createSystem(params, nonbondedMethod=5))
+        self.assertRaisesRegex(ValueError, 'LJPME is not supported', lambda:
+                    psf.createSystem(params, nonbondedMethod=LJPME))
         self.assertRaises(TypeError, lambda:
                     psf.createSystem(params, nonbondedMethod=PME,
                                      nonbondedCutoff=1*radian)
@@ -302,6 +304,9 @@ class TestCharmmFiles(unittest.TestCase):
         psf.createSystem(params, nonbondedMethod=PME, switchDistance=0.8,
                          nonbondedCutoff=1.2*nanometer)
 
+        psf_no_nbfix = CharmmPsfFile('systems/ala_ala_ala.psf', unitCellDimensions=Vec3(30, 30, 30)*angstroms)
+        psf_no_nbfix.createSystem(self.params, nonbondedMethod=LJPME)
+
     def test_ImplicitSolventForces(self):
         """Compute forces for different implicit solvent types, and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
         solventType = [HCT, OBC1, OBC2, GBn, GBn2]

wrappers/python/tests/TestDesmondDMSFile.py

@@ -38,6 +38,18 @@ class TestDesmondDMSFile(unittest.TestCase):
                                 f.getNonbondedMethod()==methodMap[method] 
                                 for f in forces))
 
+    def test_LJPME_OPLS(self):
+        """Test LJPME with and without OPLS mixing rules."""
+
+        system = self.dms_opls1.createSystem(nonbondedMethod=LJPME, OPLS=False)
+        forces = system.getForces()
+        self.assertTrue(any(isinstance(f, NonbondedForce) and
+                            f.getNonbondedMethod()==NonbondedForce.LJPME
+                            for f in forces))
+
+        with self.assertRaisesRegex(ValueError, 'LJPME is not supported'):
+            self.dms_opls1.createSystem(nonbondedMethod=LJPME, OPLS=True)
+
     def test_Cutoff(self):
         """Test to make sure the nonbondedCutoff parameter is passed correctly."""
 

wrappers/python/tests/TestForceField.py

@@ -1314,6 +1314,10 @@ class TestForceField(unittest.TestCase):
         ene2 = state2.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
         self.assertAlmostEqual(ene, ene2)
 
+        # LJPME should be forbidden with LennardJonesForce since it makes a CustomNonbondedForce to handle NBFix
+        with self.assertRaisesRegex(ValueError, 'LJPME is not supported'):
+            ff.createSystem(pdb.topology, nonbondedMethod=LJPME)
+
     def test_NBFix(self):
         """Test using LennardJonesGenerator to implement NBFix terms."""
         # Create a chain of seven atoms.

wrappers/python/tests/TestGromacsTopFile.py

@@ -36,6 +36,36 @@ class TestGromacsTopFile(unittest.TestCase):
                                 f.getNonbondedMethod()==methodMap[method]
                                 for f in forces))
 
+    def test_LJPME_mixing(self):
+        """Test that LJPME is not allowed with non-standard LJ mixing or NBFIX."""
+
+        tests = {
+           'apgr.nbfix.pairs.comb1.top': False,
+           'apgr.nbfix.pairs.comb2.top': False,
+           'apgr.nbfix.pairs.comb3.top': False,
+           'apgr.nbfix.nopairs.comb1.top': False,
+           'apgr.nbfix.nopairs.comb2.top': False,
+           'apgr.nbfix.nopairs.comb3.top': False,
+           'apgr.nonbfix.pairs.comb1.top': False,
+           'apgr.nonbfix.pairs.comb2.top': True,
+           'apgr.nonbfix.pairs.comb3.top': False,
+           'apgr.nonbfix.nopairs.comb1.top': False,
+           'apgr.nonbfix.nopairs.comb2.top': True,
+           'apgr.nonbfix.nopairs.comb3.top': False,
+        }
+
+        for topfile, ljpme_allowed in tests.items():
+            top = GromacsTopFile(f'systems/{topfile}', unitCellDimensions=Vec3(30, 30, 30)*angstroms)
+            if ljpme_allowed:
+                system = top.createSystem(nonbondedMethod=LJPME)
+                forces = system.getForces()
+                self.assertTrue(any(isinstance(f, NonbondedForce) and
+                                    f.getNonbondedMethod()==NonbondedForce.LJPME
+                                    for f in forces))
+            else:
+                with self.assertRaisesRegex(ValueError, 'LJPME is not supported'):
+                    top.createSystem(nonbondedMethod=LJPME)
+
     def test_ff99SBILDN(self):
         """ Test Gromacs topology #define replacement as used in ff99SB-ILDN """
         top = GromacsTopFile('systems/aidilnaaaaa.top')