Merge pull request #532 from swails/lj1264

Support new 12-6-4 potential

wrappers/python/CMakeLists.txt

@@ -55,6 +55,7 @@ foreach(SUBDIR ${SUBDIRS})
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdb"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop"
+        "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.inpcrd"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.parm7"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.rst7"
         "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.ncrst"

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

@@ -702,6 +702,10 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
         warnings.warn("1-4 scaling parameters in topology file are being ignored. "
             "This is not recommended unless you know what you are doing.")
 
+    has_1264 = 'LENNARD_JONES_CCOEF' in prmtop._raw_data.keys()
+    if has_1264:
+        parm_ccoef = [float(x) for x in prmtop._raw_data['LENNARD_JONES_CCOEF']]
+
     # Use pyopenmm implementation of OpenMM by default.
     if mm is None:
         mm = simtk.openmm
@@ -865,32 +869,52 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
                 idx = nbidx[numTypes*i+j] - 1
                 acoef[i+numTypes*j] = sqrt(parm_acoef[idx]) * afac
                 bcoef[i+numTypes*j] = parm_bcoef[idx] * bfac
-        cforce = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r^6;'
-                                         'a=acoef(type1, type2);'
-                                         'b=bcoef(type1, type2);')
+        if has_1264:
+            cfac = ene_conv * length_conv**4
+            ccoef = [0 for i in range(numTypes*numTypes)]
+            for i in range(numTypes):
+                for j in range(numTypes):
+                    idx = nbidx[numTypes*i+j] - 1
+                    ccoef[i+numTypes*j] = parm_ccoef[idx] * cfac
+            cforce = mm.CustomNonbondedForce('(a/r6)^2-b/r6-c/r^4; r6=r^6;'
+                                             'a=acoef(type1, type2);'
+                                             'b=bcoef(type1, type2);'
+                                             'c=ccoef(type1, type2);')
+        else:
+            cforce = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r^6;'
+                                             'a=acoef(type1, type2);'
+                                             'b=bcoef(type1, type2);')
         cforce.addTabulatedFunction('acoef',
                     mm.Discrete2DFunction(numTypes, numTypes, acoef))
         cforce.addTabulatedFunction('bcoef',
                     mm.Discrete2DFunction(numTypes, numTypes, bcoef))
+        if has_1264:
+            cforce.addTabulatedFunction('ccoef',
+                        mm.Discrete2DFunction(numTypes, numTypes, ccoef))
         cforce.addPerParticleParameter('type')
         for atom in prmtop._getAtomTypeIndexes():
             cforce.addParticle((atom-1,))
-        # Now set the various properties based on the NonbondedForce object
-        if nonbondedMethod in ('PME', 'Ewald', 'CutoffPeriodic'):
-            cforce.setNonbondedMethod(cforce.CutoffPeriodic)
-            cforce.setCutoffDistance(nonbondedCutoff)
-            cforce.setUseLongRangeCorrection(True)
-        elif nonbondedMethod == 'CutoffNonPeriodic':
-            cforce.setNonbondedMethod(cforce.CutoffNonPeriodic)
-            cforce.setCutoffDistance(nonbondedCutoff)
-        elif nonbondedMethod == 'NoCutoff':
-            cforce.setNonbondedMethod(cforce.NoCutoff)
-        else:
-            raise ValueError('Unrecognized cutoff option %s' % nonbondedMethod)
     else:
         for (charge, (rVdw, epsilon)) in zip(prmtop.getCharges(), nonbondTerms):
             sigma = rVdw * sigmaScale
             force.addParticle(charge, sigma, epsilon)
+        if has_1264:
+            numTypes = prmtop.getNumTypes()
+            nbidx = [int(x) for x in prmtop._raw_data['NONBONDED_PARM_INDEX']]
+            ccoef = [0 for i in range(numTypes*numTypes)]
+            ene_conv = units.kilocalories_per_mole.conversion_factor_to(units.kilojoules_per_mole)
+            length_conv = units.angstroms.conversion_factor_to(units.nanometers)
+            cfac = ene_conv * length_conv**4
+            for i in range(numTypes):
+                for j in range(numTypes):
+                    idx = nbidx[numTypes*i+j] - 1
+                    ccoef[i+numTypes*j] = parm_ccoef[idx] * cfac
+            cforce = mm.CustomNonbondedForce('-c/r^4; c=ccoef(type1, type2)')
+            cforce.addTabulatedFunction('ccoef',
+                        mm.Discrete2DFunction(numTypes, numTypes, ccoef))
+            cforce.addPerParticleParameter('type')
+            for atom in prmtop._getAtomTypeIndexes():
+                cforce.addParticle((atom-1,))
 
 
     # Add 1-4 Interactions
@@ -916,10 +940,23 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
 
     # Copy the exceptions as exclusions to the CustomNonbondedForce if we have
     # NBFIX terms
-    if nbfix:
+    if nbfix or has_1264:
         for i in range(force.getNumExceptions()):
             ii, jj, chg, sig, eps = force.getExceptionParameters(i)
             cforce.addExclusion(ii, jj)
+        # Now set the various properties based on the NonbondedForce object
+        if nonbondedMethod in ('PME', 'Ewald', 'CutoffPeriodic'):
+            cforce.setNonbondedMethod(cforce.CutoffPeriodic)
+            cforce.setCutoffDistance(nonbondedCutoff)
+            cforce.setUseLongRangeCorrection(True)
+        elif nonbondedMethod == 'CutoffNonPeriodic':
+            cforce.setNonbondedMethod(cforce.CutoffNonPeriodic)
+            cforce.setCutoffDistance(nonbondedCutoff)
+        elif nonbondedMethod == 'NoCutoff':
+            cforce.setNonbondedMethod(cforce.NoCutoff)
+        else:
+            raise ValueError('Unrecognized cutoff option %s' % nonbondedMethod)
+        # Add this force to the system
         system.addForce(cforce)
     system.addForce(force)
 

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -8,7 +8,9 @@ import simtk.openmm.app.element as elem
 prmtop1 = AmberPrmtopFile('systems/alanine-dipeptide-explicit.prmtop')
 prmtop2 = AmberPrmtopFile('systems/alanine-dipeptide-implicit.prmtop')
 prmtop3 = AmberPrmtopFile('systems/ff14ipq.parm7')
+prmtop4 = AmberPrmtopFile('systems/Mg_water.prmtop')
 inpcrd3 = AmberInpcrdFile('systems/ff14ipq.rst7')
+inpcrd4 = AmberInpcrdFile('systems/Mg_water.inpcrd')
 
 class TestAmberPrmtopFile(unittest.TestCase):
 
@@ -206,5 +208,38 @@ class TestAmberPrmtopFile(unittest.TestCase):
         # Amber using this force field.
         self.assertAlmostEqual(-7042.3903307/ene, 1, places=3)
 
+    def test_LJ1264(self):
+        """Test prmtop with 12-6-4 vdW potential implemented"""
+        system = prmtop4.createSystem(nonbondedMethod=PME,
+                                      nonbondedCutoff=8*angstroms)
+        # Check the forces
+        has_nonbond_force = has_custom_nonbond_force = False
+        nonbond_exceptions = custom_nonbond_exclusions = 0
+        for force in system.getForces():
+            if isinstance(force, NonbondedForce):
+                has_nonbond_force = True
+                nonbond_exceptions = force.getNumExceptions()
+                force.setUseDispersionCorrection(False)
+            elif isinstance(force, CustomNonbondedForce):
+                self.assertTrue(force.getUseLongRangeCorrection())
+                has_custom_nonbond_force = True
+                custom_nonbond_exceptions = force.getNumExclusions()
+                force.setUseLongRangeCorrection(False)
+        self.assertTrue(has_nonbond_force)
+        self.assertTrue(has_custom_nonbond_force)
+        self.assertEqual(nonbond_exceptions, custom_nonbond_exceptions)
+        integrator = VerletIntegrator(1.0*femtoseconds)
+        # Use reference platform, since it should always be present and
+        # 'working', and the system is plenty small so this won't be too slow
+        sim = Simulation(prmtop4.topology, system, integrator, Platform.getPlatformByName('Reference'))
+        # Check that the energy is about what we expect it to be
+        sim.context.setPeriodicBoxVectors(*inpcrd4.boxVectors)
+        sim.context.setPositions(inpcrd4.positions)
+        ene = sim.context.getState(getEnergy=True, enforcePeriodicBox=True).getPotentialEnergy()
+        ene = ene.value_in_unit(kilocalories_per_mole)
+        # Make sure the energy is relatively close to the value we get with
+        # Amber using this force field.
+        self.assertAlmostEqual(-7307.2735621/ene, 1, places=3)
+
 if __name__ == '__main__':
     unittest.main()