Fix reading of Gromacs topologies created from CHARMM force field (#5026)

* fix to correctly read Gromacs topology files for CHARMM force field

* updated fix to correctly read topology files with NBFIX and different combination rules

* fixes to read topology files with NBFIX and different combination rules

* changed default for useDispersionCorrection to True

* changed docstring default for useDispersionCorrection to 'True'

wrappers/python/openmm/app/gromacstopfile.py

@@ -357,7 +357,7 @@ class GromacsTopFile(object):
         if len(fields) < 3:
             raise ValueError('Too few fields in [ bonds ] line: '+line)
         if fields[2] not in ('1', '2'):
-                raise ValueError('Unsupported function type in [ bonds ] line: '+line)
+            raise ValueError('Unsupported function type in [ bonds ] line: '+line)
         self._currentMoleculeType.bonds.append(fields)
 
     def _processAngle(self, line):
@@ -652,7 +652,8 @@ class GromacsTopFile(object):
                     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, ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None, switchDistance=None):
+                     rigidWater=True, ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None, 
+                     switchDistance=None, useDispersionCorrection=True):
         """Construct an OpenMM System representing the topology described by this
         top file.
 
@@ -683,6 +684,8 @@ class GromacsTopFile(object):
         switchDistance : float=None
             The distance at which the potential energy switching function is turned on for
             Lennard-Jones interactions. If this is None, no switching function will be used.
+        useDispersionCorrection : bool=True
+            Long-range correction for dispersion (Lennard-Jones) interaction
 
         Returns
         -------
@@ -700,6 +703,19 @@ class GromacsTopFile(object):
                     atom_types.append(atom[1])
         has_nbfix_terms = any([pair in self._nonbondTypes for pair in combinations_with_replacement(sorted(set(atom_types)), 2)])
 
+        if has_nbfix_terms:
+            self._matchingNBFIX = {
+    		p: self._nonbondTypes[p]
+    		for p in combinations_with_replacement(sorted(set(atom_types)), 2)
+    		if p in self._nonbondTypes
+	    }
+            if self._defaults[1] == '3':
+                for v in self._matchingNBFIX.values():
+                    sigma=float(v[3])
+                    epsilon=float(v[4])
+                    v[3]=4*epsilon*sigma**6
+                    v[4]=4*epsilon*sigma**12
+
         # Create the System.
 
         sys = mm.System()
@@ -710,16 +726,21 @@ class GromacsTopFile(object):
             raise ValueError('Illegal nonbonded method for a non-periodic system')
         nb = mm.NonbondedForce()
         sys.addForce(nb)
-        lj = None
+       
+        lj = None        
+        ljnbfix = None   
         if has_nbfix_terms:
-            lj = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r^6; a=acoef(type1, type2); b=bcoef(type1, type2)')
-            lj.addPerParticleParameter('type')
-            sys.addForce(lj)
-        elif self._defaults[1] in ('1', '3'):
-            lj = mm.CustomNonbondedForce('A1*A2/r^12-C1*C2/r^6')
-            lj.addPerParticleParameter('C')
-            lj.addPerParticleParameter('A')
-            sys.addForce(lj)
+            # lookup table to implement NBFIX, used for all combination rules
+            ljnbfix = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r^6; a=acoef(type1, type2); b=bcoef(type1, type2)')
+            ljnbfix.addPerParticleParameter('type')
+            sys.addForce(ljnbfix)
+        else:
+            # used for LJ interactions when combination rules are 1 or 3
+            if self._defaults[1] in ('1', '3'):
+                lj = mm.CustomNonbondedForce('A1*A2/r^12-C1*C2/r^6')
+                lj.addPerParticleParameter('C')
+                lj.addPerParticleParameter('A')
+                sys.addForce(lj)
         bonds = {}
         angles = {}
         periodic = None
@@ -996,10 +1017,11 @@ class GromacsTopFile(object):
                     if has_nbfix_terms:
                         nb.addParticle(q, 1.0, 0.0)
                         atom_charges.append(q)
-                        lj.addParticle([0])
+                        ljnbfix.addParticle([0])
                     else:
                         if self._defaults[1] == '1':
                             nb.addParticle(q, 1.0, 0.0)
+                            # LJ interactions are handled via separate LJ force with custom potential
                             lj.addParticle([math.sqrt(float(params[6])), math.sqrt(float(params[7]))])
                         elif self._defaults[1] == '2':
                             nb.addParticle(q, float(params[6]), float(params[7]))
@@ -1007,6 +1029,7 @@ class GromacsTopFile(object):
                             nb.addParticle(q, 1.0, 0.0)
                             sigma = float(params[6])
                             epsilon = float(params[7])
+                            # LJ interactions are handled via separate LJ force with custom potential
                             lj.addParticle([math.sqrt(4*epsilon*sigma**6), math.sqrt(4*epsilon*sigma**12)])
 
                 for fields in moleculeType.bonds:
@@ -1035,7 +1058,7 @@ class GromacsTopFile(object):
                             atom_partners[baseAtomIndex+pair[1]]['torsion'].add(baseAtomIndex+pair[0])
 
                 # Record nonbonded exceptions.
-
+                # typically these are 1-4 pairs
                 for fields in moleculeType.pairs:
                     atoms = [int(x)-1 for x in fields[:2]]
                     types = tuple(atomTypes[i] for i in atoms)
@@ -1043,12 +1066,24 @@ class GromacsTopFile(object):
                     atom2params = nb.getParticleParameters(baseAtomIndex+atoms[1])
                     atom1params = [x.value_in_unit_system(unit.md_unit_system) for x in atom1params]
                     atom2params = [x.value_in_unit_system(unit.md_unit_system) for x in atom2params]
+
+                    def convertParams(params):
+                        if self._defaults[1] == '3':
+                           # convert from sigma/epsilon given in topology file for combination rule 3 according to GROMACS convention
+                           sigma=params[0]
+                           epsilon=params[1]
+                           return [4*epsilon*sigma**6, 4*epsilon*sigma**12]
+                        return params
+                        
                     if len(fields) >= 5:
-                        params = [float(x) for x in fields[3:5]]
+                        # extra parameters given for 1-4 interactions as part of the pair entry
+                        params = convertParams([float(x) for x in fields[3:5]])
                     elif types in self._pairTypes:
-                        params = [float(x) for x in self._pairTypes[types][3:5]]
+                        # parameters given under pairTypes
+                        params = convertParams([float(x) for x in self._pairTypes[types][3:5]])
                     elif types[::-1] in self._pairTypes:
-                        params = [float(x) for x in self._pairTypes[types[::-1]][3:5]]
+                        # parameters given under pairTypes
+                        params = convertParams([float(x) for x in self._pairTypes[types[::-1]][3:5]])
                     elif not self._genpairs:
                         raise ValueError('No pair parameters defined for atom '
                                          'types %s and gen-pairs is "no"' % types)
@@ -1125,10 +1160,15 @@ class GromacsTopFile(object):
                 charge_prod = fudgeQQ*q1*q4
                 epsilon = math.sqrt(abs(eps1 * eps4))
                 if self._defaults[1] == '2':
-                   rmin14 = (rmin1 + rmin4) / 2
+                    rmin14 = (rmin1 + rmin4) / 2
                 else:
-                   rmin14 = math.sqrt(rmin1 * rmin4)
-                nb.addException(tor[0], tor[3], charge_prod, rmin14, epsilon)
+                    rmin14 = math.sqrt(rmin1 * rmin4)
+                # Parameters are generated via standard combining rules.
+                # If different 1-4 parameters are given via pairtypes they will be overwritten below.
+                if self._defaults[1] == '3':
+                    nb.addException(tor[0], tor[3], charge_prod, 4*epsilon*rmin14**6, 4*epsilon*rmin14**12)
+                else:
+                    nb.addException(tor[0], tor[3], charge_prod, rmin14, epsilon)
                 excluded_atom_pairs.add(key)
 
             # Add excluded atoms
@@ -1153,25 +1193,33 @@ class GromacsTopFile(object):
         for exclusion in exclusions:
             nb.addException(exclusion[0], exclusion[1], 0.0, 1.0, 0.0, True)
 
-        if lj is not None:
-            # We're using a CustomNonbondedForce for LJ interactions, so also create a CustomBondForce
-            # to handle the exceptions.
+        # this will overwrite the pairs from the pairlist
+        # if nbfix, this will only overwrite paris if we have pairtype parameters 
+        for pair in pairs:
+            nb.addException(pair[0], pair[1], pair[2], pair[3], pair[4], True)
 
+        if self._defaults[1] in ('1', '3'):
+           # We're using a CustomNonbondedForce for LJ interactions, so also create a CustomBondForce
+           # to handle the exceptions
+ 
             pair_bond = mm.CustomBondForce('-C/r^6+A/r^12')
             pair_bond.addPerBondParameter('C')
             pair_bond.addPerBondParameter('A')
             pair_bond.setName('LennardJonesExceptions')
             sys.addForce(pair_bond)
-            for pair in pairs:
-                nb.addException(pair[0], pair[1], pair[2], 1.0, 0.0, True)
-                pair_bond.addBond(pair[0], pair[1], [pair[3], pair[4]])
+            for i in range(nb.getNumExceptions()):
+                ii, jj, q, sig, eps = nb.getExceptionParameters(i)
+                if (eps.value_in_unit(unit.kilojoule_per_mole))>0:
+                    nb.addException(ii,jj,q, 1.0, 0.0, True)
+                    pair_bond.addBond(ii,jj,[sig,eps])
+                if lj is not None:
+                    lj.addExclusion(ii, jj)
+
+        if ljnbfix is not None:
             for i in range(nb.getNumExceptions()):
                 ii, jj, q, eps, sig = nb.getExceptionParameters(i)
-                lj.addExclusion(ii, jj)
-        elif self._defaults[1] == '2':
-            for pair in pairs:
-                nb.addException(pair[0], pair[1], pair[2], pair[3], pair[4], True)
-
+                ljnbfix.addExclusion(ii, jj)
+           
         # Finish configuring the NonbondedForce.
 
         methodMap = {ff.NoCutoff:mm.NonbondedForce.NoCutoff,
@@ -1183,9 +1231,14 @@ class GromacsTopFile(object):
         nb.setNonbondedMethod(methodMap[nonbondedMethod])
         nb.setCutoffDistance(nonbondedCutoff)
         nb.setEwaldErrorTolerance(ewaldErrorTolerance)
+
+        if nonbondedMethod in (ff.PME, ff.LJPME, ff.Ewald, ff.CutoffPeriodic):
+            nb.setUseDispersionCorrection(bool(useDispersionCorrection))
+
         if switchDistance is not None:
             nb.setUseSwitchingFunction(True)
             nb.setSwitchingDistance(switchDistance)
+
         if lj is not None:
             methodMap = {ff.NoCutoff:mm.CustomNonbondedForce.NoCutoff,
                          ff.CutoffNonPeriodic:mm.CustomNonbondedForce.CutoffNonPeriodic,
@@ -1195,16 +1248,29 @@ class GromacsTopFile(object):
                          ff.LJPME:mm.CustomNonbondedForce.CutoffPeriodic}
             lj.setNonbondedMethod(methodMap[nonbondedMethod])
             lj.setCutoffDistance(nonbondedCutoff)
-            if nonbondedMethod in (ff.PME, ff.LJPME, ff.Ewald, ff.CutoffPeriodic):
-                lj.setUseLongRangeCorrection(True)
+            lj.setUseLongRangeCorrection(bool(useDispersionCorrection))
             if switchDistance is not None:
                 lj.setUseSwitchingFunction(True)
                 lj.setSwitchingDistance(switchDistance)
             lj.setName('LennardJonesForce')
 
         if has_nbfix_terms:
+            methodMap = {ff.NoCutoff:mm.CustomNonbondedForce.NoCutoff,
+                         ff.CutoffNonPeriodic:mm.CustomNonbondedForce.CutoffNonPeriodic,
+                         ff.CutoffPeriodic:mm.CustomNonbondedForce.CutoffPeriodic,
+                         ff.Ewald:mm.CustomNonbondedForce.CutoffPeriodic,
+                         ff.PME:mm.CustomNonbondedForce.CutoffPeriodic,
+                         ff.LJPME:mm.CustomNonbondedForce.CutoffPeriodic}
+            ljnbfix.setNonbondedMethod(methodMap[nonbondedMethod])
+            ljnbfix.setCutoffDistance(nonbondedCutoff)
+            ljnbfix.setUseLongRangeCorrection(bool(useDispersionCorrection))
+            if switchDistance is not None:
+                ljnbfix.setUseSwitchingFunction(True)
+                ljnbfix.setSwitchingDistance(switchDistance)
+            ljnbfix.setName('LennardJonesForce')
+
             atom_nbfix_types = set([])
-            for pair in self._nonbondTypes:
+            for pair in self._matchingNBFIX:
                 atom_nbfix_types.add(pair[0])
                 atom_nbfix_types.add(pair[1])
 
@@ -1228,7 +1294,7 @@ class GromacsTopFile(object):
                     ljtype2 = (float(atom2[6]), float(atom2[7]))
                     if atom2 is atom:
                         lj_idx_list[j] = num_lj_types
-                    elif atom_type not in atom_nbfix_types:
+                    elif (atom_type not in atom_nbfix_types) and (atom_type2 not in atom_nbfix_types):
                         # Only non-NBFIXed atom types can be compressed
                         if ljtype == ljtype2:
                             lj_idx_list[j] = num_lj_types
@@ -1242,7 +1308,7 @@ class GromacsTopFile(object):
                 for j in range(num_lj_types):
                     namej = lj_type_list[j][0]
                     try:
-                        types = self._nonbondTypes[tuple(sorted((namei, namej)))]
+                        types = self._matchingNBFIX[tuple(sorted((namei, namej)))]
                         params = (float(types[3]), float(types[4]))
                         if self._defaults[1] == '2':
                             c6 = 4 * params[1] * params[0]**6
@@ -1260,16 +1326,16 @@ class GromacsTopFile(object):
                             if self._defaults[1] == '2':
                                 sigma = (params1[0] + params2[0]) / 2
                             else:
-                                sigma = math.sqrt(params1[0] + params2[0])
+                                sigma = math.sqrt(params1[0] * params2[0])
                             epsilon = math.sqrt(params1[1] * params2[1])
                             c6 = 4 * epsilon * sigma**6
                             c12 = 4 * epsilon * sigma**12
                     acoef[i+num_lj_types*j] = math.sqrt(c12)
                     bcoef[i+num_lj_types*j] = c6
-            lj.addTabulatedFunction('acoef', mm.Discrete2DFunction(num_lj_types, num_lj_types, acoef))
-            lj.addTabulatedFunction('bcoef', mm.Discrete2DFunction(num_lj_types, num_lj_types, bcoef))
+            ljnbfix.addTabulatedFunction('acoef', mm.Discrete2DFunction(num_lj_types, num_lj_types, acoef))
+            ljnbfix.addTabulatedFunction('bcoef', mm.Discrete2DFunction(num_lj_types, num_lj_types, bcoef))
             for i, idx in enumerate(lj_idx_list):
-                lj.setParticleParameters(i, [idx-1]) # adjust for indexing from 0
+                ljnbfix.setParticleParameters(i, [idx-1]) # adjust for indexing from 0
 
         # Adjust masses.
 

wrappers/python/tests/TestGromacsTopFile.py

@@ -213,7 +213,7 @@ class TestGromacsTopFile(unittest.TestCase):
             if i == 4:
                 wc = -wc
             self.assertAlmostEqual(wc, vs.getWeightCross())
-
+    
     def test_GROMOS(self):
         """Test a system using the GROMOS 54a7 force field."""
 
@@ -240,6 +240,34 @@ class TestGromacsTopFile(unittest.TestCase):
         total = context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilojoules_per_mole)
         assert_allclose(-1.77020e+02, total, rtol=1e-3)
 
+    def test_NBFIX(self):
+        """Test systems using NBFIX with and without pairtypes and different combination rules."""
+
+        gro = GromacsGroFile('systems/apgr.gro')
+
+        energies = {
+           'apgr.nbfix.pairs.comb1.top': -986.713,
+           'apgr.nbfix.pairs.comb2.top': -986.485,
+           'apgr.nbfix.pairs.comb3.top': -986.713,
+           'apgr.nbfix.nopairs.comb1.top': -912.181,
+           'apgr.nbfix.nopairs.comb2.top': -903.437,
+           'apgr.nbfix.nopairs.comb3.top': -912.181,
+           'apgr.nonbfix.pairs.comb1.top': -993.104,
+           'apgr.nonbfix.pairs.comb2.top': -992.685,
+           'apgr.nonbfix.pairs.comb3.top': -993.104,
+           'apgr.nonbfix.nopairs.comb1.top': -918.572,
+           'apgr.nonbfix.nopairs.comb2.top': -909.637,
+           'apgr.nonbfix.nopairs.comb3.top': -918.572
+        }
+        
+        for topfile, expected in energies.items():   
+           top = GromacsTopFile(f'systems/{topfile}')
+           system = top.createSystem(nonbondedMethod=NoCutoff,switchDistance=None,constraints=None,useDispersionCorrection=False)
+           context = Context(system, VerletIntegrator(1*femtosecond), Platform.getPlatform('Reference'))
+           context.setPositions(gro.positions)
+           energy=context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilojoules_per_mole)
+           assert_allclose(energy, expected, rtol=1E-6)
+
 if __name__ == '__main__':
     unittest.main()
 

wrappers/python/tests/systems/apgr.gro

+GROMACS
+   58
+    1ALA      N    1   2.374   2.436   2.592
+    1ALA     H1    2   2.456   2.469   2.535
+    1ALA     H2    3   2.292   2.473   2.541
+    1ALA     H3    4   2.375   2.332   2.583
+    1ALA     CA    5   2.389   2.483   2.732
+    1ALA     HA    6   2.431   2.583   2.729
+    1ALA     CB    7   2.249   2.486   2.798
+    1ALA    HB1    8   2.204   2.384   2.802
+    1ALA    HB2    9   2.180   2.553   2.742
+    1ALA    HB3   10   2.256   2.525   2.902
+    1ALA      C   11   2.485   2.395   2.810
+    1ALA      O   12   2.505   2.280   2.774
+    2PRO      N   13   2.545   2.453   2.915
+    2PRO     CD   14   2.525   2.378   3.041
+    2PRO    HD1   15   2.439   2.307   3.038
+    2PRO    HD2   16   2.505   2.453   3.121
+    2PRO     CA   17   2.690   2.474   2.909
+    2PRO     HA   18   2.714   2.547   2.985
+    2PRO     CB   19   2.751   2.341   2.955
+    2PRO    HB1   20   2.856   2.352   2.989
+    2PRO    HB2   21   2.744   2.263   2.876
+    2PRO     CG   22   2.657   2.304   3.071
+    2PRO    HG1   23   2.642   2.194   3.077
+    2PRO    HG2   24   2.699   2.340   3.168
+    2PRO      C   25   2.747   2.535   2.782
+    2PRO      O   26   2.759   2.657   2.775
+    3GLY      N   27   2.778   2.453   2.679
+    3GLY     HN   28   2.770   2.355   2.693
+    3GLY     CA   29   2.809   2.499   2.544
+    3GLY    HA1   30   2.856   2.417   2.492
+    3GLY    HA2   31   2.872   2.587   2.551
+    3GLY      C   32   2.685   2.539   2.467
+    3GLY      O   33   2.589   2.591   2.523
+    4ARG      N   34   2.681   2.511   2.336
+    4ARG     HN   35   2.758   2.465   2.292
+    4ARG     CA   36   2.566   2.535   2.250
+    4ARG     HA   37   2.553   2.643   2.245
+    4ARG     CB   38   2.597   2.481   2.108
+    4ARG    HB1   39   2.623   2.372   2.119
+    4ARG    HB2   40   2.689   2.532   2.071
+    4ARG     CG   41   2.487   2.491   1.999
+    4ARG    HG1   42   2.402   2.424   2.027
+    4ARG    HG2   43   2.531   2.448   1.906
+    4ARG     CD   44   2.435   2.631   1.965
+    4ARG    HD1   45   2.380   2.628   1.868
+    4ARG    HD2   46   2.518   2.704   1.958
+    4ARG     NE   47   2.343   2.670   2.076
+    4ARG     HE   48   2.361   2.627   2.166
+    4ARG     CZ   49   2.215   2.702   2.065
+    4ARG    NH1   50   2.141   2.679   2.171
+    4ARG   HH11   51   2.189   2.630   2.245
+    4ARG   HH12   52   2.043   2.688   2.167
+    4ARG    NH2   53   2.163   2.754   1.957
+    4ARG   HH21   54   2.068   2.784   1.959
+    4ARG   HH22   55   2.226   2.781   1.885
+    4ARG      C   56   2.434   2.474   2.305
+    4ARG    OT1   57   2.431   2.353   2.336
+    4ARG    OT2   58   2.336   2.551   2.326
+   5.00000   5.00000   5.00000