Merge pull request #2120 from jing-huang/drude

Implementing Drude polarizable force fields in the CHARMM file parsers

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

@@ -12,7 +12,7 @@ the ParmEd program and was ported for use with OpenMM.
 Copyright (c) 2014 the Authors
 
 Author: Jason M. Swails
-Contributors:
+Contributors: Jing Huang
 Date: Sep. 17, 2014
 
 Permission is hereby granted, free of charge, to any person obtaining a
@@ -114,6 +114,7 @@ class CharmmParameterSet(object):
         self.improper_types = dict()
         self.cmap_types = dict()
         self.nbfix_types = dict()
+        self.nbthole_types = dict()
         self.parametersets = []
 
         # Load all of the files
@@ -261,6 +262,9 @@ class CharmmParameterSet(object):
             if line.startswith('NBFIX'):
                 section = 'NBFIX'
                 continue
+            if line.startswith('THOLE'):
+                section = 'NBTHOLE'
+                continue
             if line.startswith('HBOND'):
                 section = None
                 continue
@@ -493,6 +497,22 @@ class CharmmParameterSet(object):
                 except IndexError:
                     raise CharmmFileError('Could not parse NBFIX terms.')
                 self.nbfix_types[(min(at1, at2), max(at1, at2))] = (emin, rmin)
+                continue
+            # Here parse the possible nbthole section
+            if section == 'NBTHOLE':
+                words = line.split()
+                try:
+                    at1 = words[0]
+                    at2 = words[1]
+                    nbt = abs(conv(words[2], float, 'NBTHOLE a'))
+                    try:
+                        self.atom_types_str[at1].add_nbthole(at2, nbt)
+                        self.atom_types_str[at2].add_nbthole(at1, nbt)
+                    except KeyError:
+                        pass
+                except IndexError:
+                    raise CharmmFileError('Could not parse NBTHOLE terms.')
+                self.nbthole_types[(min(at1, at2), max(at1, at2))] = (nbt)
         # If there were any CMAP terms stored in the parameter set, the last one
         # defined will not have been added to the set. Add it now.
         if current_cmap is not None:

wrappers/python/simtk/openmm/app/internal/charmm/topologyobjects.py

@@ -100,6 +100,9 @@ class AtomType(object):
     nbfix : dict
         Dictionary that maps nbfix terms with other atom types. Dict entries
         are (rmin, epsilon) -- precombined values for that particular atom pair
+    nbthole : dict
+            Dictionary that maps nbthole terms with other atom types. Dict entries
+            are the value of a -- screening factor for that pair
 
     Examples
     --------
@@ -137,6 +140,8 @@ class AtomType(object):
         # Store each NBFIX term as a dict with the atom type string matching to
         # a 2-element tuple that is rmin, epsilon
         self.nbfix = dict()
+        # Likewise, store each NBTHOLE term as a dict 
+        self.nbthole = dict()
 
     def __eq__(self, other):
         """
@@ -173,6 +178,10 @@ class AtomType(object):
         if epsilon14 is None: epsilon14 = epsilon
         self.nbfix[typename] = (rmin, epsilon, rmin14, epsilon14)
 
+    def add_nbthole(self, typename, nbt):
+        """ Adds a new NBTHOLE screening factor for this atom """
+        self.nbthole[typename] = (nbt)
+
     def __str__(self):
         return self.name
 

wrappers/python/tests/TestCharmmFiles.py

@@ -116,6 +116,26 @@ class TestCharmmFiles(unittest.TestCase):
         ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
         self.assertAlmostEqual(ene, 15490.0033559, delta=0.05)
 
+    def test_drude(self):
+        """Tests CHARMM systems with Drude force field"""
+        warnings.filterwarnings('ignore', category=CharmmPSFWarning)
+        psf = CharmmPsfFile('systems/ala3_solv_drude.psf')
+        crd = CharmmCrdFile('systems/ala3_solv_drude.crd')
+        params = CharmmParameterSet('systems/toppar_drude_master_protein_2013e.str')
+        # Box dimensions (cubic box)
+        psf.setBox(33.2*angstroms, 33.2*angstroms, 33.2*angstroms)
+
+        # Now compute the full energy
+        plat = Platform.getPlatformByName('Reference')
+        system = psf.createSystem(params, nonbondedMethod=PME)
+        integrator = DrudeLangevinIntegrator(300*kelvin, 1.0/picosecond, 1*kelvin, 10/picosecond, 0.001*picoseconds)
+        con = Context(system, integrator, plat)
+        con.setPositions(crd.positions)
+
+        state = con.getState(getEnergy=True, enforcePeriodicBox=True)
+        ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
+        self.assertAlmostEqual(ene, -1831.54, delta=0.5)
+
     def test_InsCode(self):
         """ Test the parsing of PSF files that contain insertion codes in their residue numbers """
         psf = CharmmPsfFile('systems/4TVP-dmj_wat-ion.psf')