Merge pull request #1 from peastman/gromacstop

Fixes to GromacsTopFile with different combining rules

wrappers/python/simtk/openmm/app/gromacstopfile.py

@@ -596,15 +596,13 @@ class GromacsTopFile(object):
             sys.setDefaultPeriodicBoxVectors(*boxVectors)
         elif nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME, ff.LJPME):
             raise ValueError('Illegal nonbonded method for a non-periodic system')
-        if self._defaults[1] == '1':
-            # factor of 138.935456 for coulomb taken from tests/TestCustomNonbondedForce.h:574
-            nb = mm.CustomNonbondedForce('138.935456*q1*q2/r-sqrt(C1*C2)/r^6+sqrt(A1*A2)/r^12')
-            nb.addPerParticleParameter('q')
-            nb.addPerParticleParameter('C')
-            nb.addPerParticleParameter('A')
-        elif self._defaults[1] == '2':
-            nb = mm.NonbondedForce()
+        nb = mm.NonbondedForce()
         sys.addForce(nb)
+        if self._defaults[1] == '1':
+            lj = mm.CustomNonbondedForce('sqrt(A1*A2)/r^12-sqrt(C1*C2)/r^6')
+            lj.addPerParticleParameter('C')
+            lj.addPerParticleParameter('A')
+            sys.addForce(lj)
         if implicitSolvent is OBC2:
             gb = mm.GBSAOBCForce()
             gb.setSoluteDielectric(soluteDielectric)
@@ -622,7 +620,7 @@ class GromacsTopFile(object):
         mapIndices = {}
         bondIndices = []
         topologyAtoms = list(self.topology.atoms())
-        exceptions = []
+        exclusions = []
         pairs = []
         fudgeQQ = float(self._defaults[4])
         fudgeLJ = float(self._defaults[3])
@@ -793,7 +791,7 @@ class GromacsTopFile(object):
                             phi0 = float(params[5])
                             if k != 0:
                                 if harmonicTorsion is None:
-                                    harmonicTorsion = mm.CustomTorsionForce('0.5*k*(thetap-theta0)^2; thetap = step(-(theta-theta0+pi))*2*pi+theta+step(theta-theta0-pi)*(-2*pi); pi = 3.14159')
+                                    harmonicTorsion = mm.CustomTorsionForce('0.5*k*(thetap-theta0)^2; thetap = step(-(theta-theta0+pi))*2*pi+theta+step(theta-theta0-pi)*(-2*pi); pi = %.15g' % math.pi)
                                     harmonicTorsion.addPerTorsionParameter('theta0')
                                     harmonicTorsion.addPerTorsionParameter('k')
                                     sys.addForce(harmonicTorsion)
@@ -849,7 +847,8 @@ class GromacsTopFile(object):
                         q = float(params[4])
 
                     if self._defaults[1] == '1':
-                        nb.addParticle([q, float(params[6]), float(params[7])])
+                        nb.addParticle(q, 1.0, 0.0)
+                        lj.addParticle([float(params[6]), float(params[7])])
                     elif self._defaults[1] == '2':
                         nb.addParticle(q, float(params[6]), float(params[7]))
 
@@ -880,44 +879,37 @@ class GromacsTopFile(object):
                         continue # We'll use the automatically generated parameters
                     atom1params = nb.getParticleParameters(baseAtomIndex+atoms[0])
                     atom2params = nb.getParticleParameters(baseAtomIndex+atoms[1])
-                    # enable tracking pairs and exceptions separately
-                    if self._defaults[1] == '1':
-                        pairs.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1]))
-                    elif self._defaults[1] == '2':
-                        exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1]))
+                    pairs.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1]))
                 for fields in moleculeType.exclusions:
                     atoms = [int(x)-1 for x in fields]
                     for atom in atoms[1:]:
                         if atom > atoms[0]:
-                            if self._defaults[1] == '1':
-                                exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atom))
-                            elif self._defaults[1] == '2':
-                                exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atom, 0, 0, 0))
+                            exclusions.append((baseAtomIndex+atoms[0], baseAtomIndex+atom))
 
 
         # Create nonbonded exceptions.
+        
+        nb.createExceptionsFromBonds(bondIndices, fudgeQQ, fudgeLJ)
+        for exclusion in exclusions:
+            nb.addException(exclusion[0], exclusion[1], 0.0, 1.0, 0.0, True)
         if self._defaults[1] == '1':
-            # implement named exceptions (generally used for 1-4 interactions)
+            # We're using a CustomNonbondedForce for LJ interactions, so also create a CustomBondForce
+            # to handle the exceptions.
+            
             pair_bond = mm.CustomBondForce('138.935456*q/r-C/r^6+A/r^12')
             pair_bond.addPerBondParameter('q')
             pair_bond.addPerBondParameter('C')
             pair_bond.addPerBondParameter('A')
             sys.addForce(pair_bond)
-
-            # generate 1-2, 1-3, and 1-4 exclusions
-            nb.createExclusionsFromBonds(bondIndices, 3)
-
+            lj.createExclusionsFromBonds(bondIndices, 3)
             for pair in pairs:
-                pair_bond.addBond(pair[0], pair[1], [pair[2],float(pair[3])*fudgeLJ, float(pair[4])*fudgeLJ])
-
-            # create requested exclusions
-            for exception in exceptions:
-                nb.addExclusion(exception[0], exception[1])
-
+                nb.addException(pair[0], pair[1], pair[2], 1.0, 0.0, True)
+                pair_bond.addBond(pair[0], pair[1], [pair[2],float(pair[3]), float(pair[4])])
+            for exclusion in exclusions:
+                lj.addExclusion(exclusion[0], exclusion[1])
         elif self._defaults[1] == '2':
-            nb.createExceptionsFromBonds(bondIndices, fudgeQQ, fudgeLJ)
-            for exception in exceptions:
-                nb.addException(exception[0], exception[1], exception[2], float(exception[3])*fudgeLJ, float(exception[4]), True)
+            for pair in pairs:
+                nb.addException(pair[0], pair[1], pair[2], float(pair[3]), float(pair[4]), True)
 
         # Finish configuring the NonbondedForce.
 
@@ -929,8 +921,16 @@ class GromacsTopFile(object):
                      ff.LJPME:mm.NonbondedForce.LJPME}
         nb.setNonbondedMethod(methodMap[nonbondedMethod])
         nb.setCutoffDistance(nonbondedCutoff)
-        if self._defaults[1] == '2':
-            nb.setEwaldErrorTolerance(ewaldErrorTolerance)
+        nb.setEwaldErrorTolerance(ewaldErrorTolerance)
+        if self._defaults[1] == '1':
+            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}
+            lj.setNonbondedMethod(methodMap[nonbondedMethod])
+            lj.setCutoffDistance(nonbondedCutoff)
 
         # Adjust masses.