Support of Amber FF19SB CMAP (#3143)

* Support of Amber FF19SB CMAP

* Update: Support of Amber FF19SB CMAP
Co-authored-by: Kye Won Wang <kyw220@sol.cc.lehigh.edu>

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

@@ -14,7 +14,7 @@ Medical Research, grant U54 GM072970. See https://simtk.org.
 
 Portions copyright (c) 2012-2014 Stanford University and the Authors.
 Authors: Randall J. Radmer, John D. Chodera, Peter Eastman
-Contributors: Christoph Klein, Michael R. Shirts, Jason Swails
+Contributors: Christoph Klein, Michael R. Shirts, Jason Swails, Kye Won Wang
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -62,7 +62,7 @@ from . import customgbforces as customgb
 #=============================================================================================
 
 # A regex for extracting print format info from the FORMAT lines.
-FORMAT_RE_PATTERN=re.compile("([0-9]+)([a-zA-Z]+)([0-9]+)\.?([0-9]*)")
+FORMAT_RE_PATTERN=re.compile("([0-9]+)\(?([a-zA-Z]+)([0-9]+)\.?([0-9]*)\)?")
 
 # Pointer labels which map to pointer numbers at top of prmtop files
 POINTER_LABELS  = """
@@ -132,8 +132,6 @@ class PrmtopLoader(object):
                             raise TypeError('CHAMBER-style topology files are not supported here. '
                                             'Consider using the CHARMM files directly with CharmmPsfFile '
                                             'or ParmEd (where CHAMBER topologies are supported)')
-                        if 'CMAP' in flag:
-                            raise TypeError("CMAP terms in AMBER topology files are not supported. You can use ParmEd instead")
                         self._flags.append(flag)
                         self._raw_data[flag] = []
                     elif line.startswith('%FORMAT'):
@@ -461,6 +459,62 @@ class PrmtopLoader(object):
                                 int(0.5+float(periodicity[iType]))))
         return self._dihedralList
 
+    def getNumMaps(self):
+        """Return number of CMAPs. Return 0 if CMAP does not exist"""
+        try:
+            return self._numCMAP
+        except AttributeError:
+            pass
+        if "CMAP_COUNT" in self._raw_data.keys():
+            self._numCMAP=int(self._raw_data["CMAP_COUNT"][1])
+            return self._numCMAP
+        return 0
+
+    def getCMAPResolutions(self):
+        """Return CMAP resolution info. Return 0 if CMAP does not exist"""
+        try:
+            return self._cmapResolution
+        except AttributeError:
+            pass
+        if "CMAP_RESOLUTION" in self._raw_data.keys():
+            self._cmapResolution=self._raw_data["CMAP_RESOLUTION"]
+            return self._cmapResolution
+        return 0
+
+    def getCMAPParameters(self, index):
+        """Return list of CMAP energy values"""
+        flag="CMAP_PARAMETER_{:02d}".format(index)
+        return [float(pointer) for pointer in self._raw_data[flag]]
+
+    def getCMAPDihedrals(self):
+        """Return CMAP type, list of first four atoms, and list of second four atoms"""
+        try:
+            return self._cmapList
+        except AttributeError:
+            pass
+        cmapPointers = self._raw_data["CMAP_INDEX"]
+        self._cmapList=[]
+        forceConstConversionFactor = (units.kilocalorie_per_mole).conversion_factor_to(units.kilojoule_per_mole)
+        for ii in range(0,len(cmapPointers),6):
+            if any([int(cmapPointers[ii+jj])<0 for jj in range(5)]):
+                raise ValueError("Found negative cmap atom pointers %s"
+                                % ((cmapPointers[ii],
+                                   cmapPointers[ii+1],
+                                   cmapPointers[ii+2],
+                                   cmapPointers[ii+3],
+                                   cmapPointers[ii+4]),))
+            iType=int(cmapPointers[ii+5])-1
+            self._cmapList.append((int(iType),
+                               int(cmapPointers[ii])-1,
+                               int(cmapPointers[ii+1])-1,
+                               int(cmapPointers[ii+2])-1,
+                               int(cmapPointers[ii+3])-1,
+                               int(cmapPointers[ii+1])-1,
+                               int(cmapPointers[ii+2])-1,
+                               int(cmapPointers[ii+3])-1,
+                               int(cmapPointers[ii+4])-1))
+        return self._cmapList
+
     def get14Interactions(self):
         """Return list of atom pairs, chargeProduct, rMin and epsilon for each 1-4 interaction"""
         dihedralPointers = self._raw_data["DIHEDRALS_INC_HYDROGEN"] \
@@ -742,6 +796,34 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
         force.addTorsion(iAtom, jAtom, kAtom, lAtom, periodicity, phase, forceConstant)
     system.addForce(force)
 
+    # Add CMAP info.
+    ## Get mapSize and Resolutions
+    numMap = prmtop.getNumMaps()
+    mapSize = prmtop.getCMAPResolutions()
+
+    if numMap > 0:
+        if verbose: print("Adding CMAPs...")
+        force = mm.CMAPTorsionForce()
+        ### Get map energies
+        for field in range(numMap):
+            index = field + 1
+            ngrid = int(mapSize[field])
+            cmap = []
+            cmap_param = prmtop.getCMAPParameters(index)
+            forceConstConversionFactor = (units.kilocalorie_per_mole).conversion_factor_to(units.kilojoule_per_mole)
+            for i in range(ngrid):
+                for j in range(ngrid):
+                    idx = ngrid*((j+ngrid//2)%ngrid)+((i+ngrid//2)%ngrid)
+                    cmap.append(cmap_param[idx]*forceConstConversionFactor)
+            cmap = tuple(cmap)
+            force.addMap(ngrid, cmap)
+
+        #### Add CMAPtorsions.
+        if verbose: print("Adding CMAP torsions...")
+        for (Type, iAtom, jAtom, kAtom, lAtom, jAtom, kAtom, lAtom, mAtom) in prmtop.getCMAPDihedrals():
+            index=force.addTorsion(Type, iAtom, jAtom, kAtom, lAtom, jAtom, kAtom, lAtom, mAtom)
+        system.addForce(force)
+
     # Add nonbonded interactions.
     if verbose: print("Adding nonbonded interactions...")
     force = mm.NonbondedForce()

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -13,8 +13,10 @@ prmtop3 = AmberPrmtopFile('systems/ff14ipq.parm7')
 prmtop4 = AmberPrmtopFile('systems/Mg_water.prmtop')
 prmtop5 = AmberPrmtopFile('systems/tz2.truncoct.parm7')
 prmtop6 = AmberPrmtopFile('systems/gaffwat.parm7')
+prmtop7 = AmberPrmtopFile('systems/18protein.parm7')
 inpcrd3 = AmberInpcrdFile('systems/ff14ipq.rst7')
 inpcrd4 = AmberInpcrdFile('systems/Mg_water.inpcrd')
+inpcrd7 = AmberInpcrdFile('systems/18protein.rst7')
 
 class TestAmberPrmtopFile(unittest.TestCase):
 
@@ -403,5 +405,32 @@ class TestAmberPrmtopFile(unittest.TestCase):
             energy = context.getState(getEnergy=True, groups={1}).getPotentialEnergy().value_in_unit(kilojoules_per_mole)
             self.assertAlmostEqual(energy, expectedEnergy, delta=5e-4*abs(energy))
 
+    def testAmberCMAP(self):
+        """Check that CMAP energy calcultion compared to AMber."""
+        temperature = 50*kelvin
+        conversion = 4.184 # 4.184 kJ/mol
+        sander_CMAP_E = 8.2864 # CMAP energy calcluated by Amber, unit kcal/mol
+
+        prmtop = prmtop7  # systems/18protein.parm7
+        inpcrd = inpcrd7  # systems/18protein.rst7
+
+        system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1.2)
+        integrator = LangevinIntegrator(temperature, 1.0 / picosecond, 0.002 * picoseconds)
+
+        simulation = Simulation(prmtop.topology, system, integrator)
+        simulation.context.setPositions(inpcrd.positions)
+        simulation.context.setPeriodicBoxVectors(*inpcrd.boxVectors)
+        
+        for i, force in enumerate(system.getForces()):
+            force.setForceGroup(i)
+            
+        simulation.context.reinitialize(True)
+
+        for i in range(system.getNumForces()):
+            if i == 3: # 3 indicates CMAP force
+#                print(simulation.context.getState(getEnergy=True, groups=1<<i).getPotentialEnergy().value_in_unit(kilojoules_per_mole))
+                OpenMM_CMAP_E = simulation.context.getState(getEnergy=True, groups=1<<i).getPotentialEnergy().value_in_unit(kilojoules_per_mole)/conversion
+                self.assertAlmostEqual(OpenMM_CMAP_E, sander_CMAP_E, places=4)
+
 if __name__ == '__main__':
     unittest.main()