Merge pull request #2397 from peastman/gbn2

Use correct parameters for nucleic acids with GBn2

wrappers/python/simtk/openmm/app/data/pdbNames.xml

@@ -239,19 +239,19 @@
   <Atom name="HN2" alt1="2H" alt2="HT2" alt3="H2"/>
  </Residue>
  <Residue name="UNK" type="Protein"/>
- <Residue name="A" alt1="ADE" type="Nucleic">
+ <Residue name="A" alt1="ADE" alt2="A3" alt3="A5" type="Nucleic">
   <Atom name="H61" alt1="1H6"/>
   <Atom name="H62" alt1="2H6"/>
  </Residue>
- <Residue name="G" alt1="GUA" type="Nucleic">
+ <Residue name="G" alt1="GUA" alt2="G3" alt3="G5" type="Nucleic">
   <Atom name="H21" alt1="1H2"/>
   <Atom name="H22" alt1="2H2"/>
  </Residue>
- <Residue name="C" alt1="CYT" type="Nucleic">
+ <Residue name="C" alt1="CYT" alt2="C3" alt3="C5" type="Nucleic">
   <Atom name="H41" alt1="1H4"/>
   <Atom name="H42" alt1="2H4"/>
  </Residue>
- <Residue name="U" alt1="URA" type="Nucleic"/>
+ <Residue name="U" alt1="URA" alt2="U3" alt3="U5" type="Nucleic"/>
  <Residue name="DA" alt1="DAD" alt2="DA3" alt3="DA5" type="Nucleic">
   <Atom name="H61" alt1="1H6"/>
   <Atom name="H62" alt1="2H6"/>

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

@@ -389,18 +389,20 @@ def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, k
                                 "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))"+params, CustomGBForce.ParticlePairNoExclusions)
 
 
-_SCREEN_PARAMETERS = { # normal, GBn, GBn2
-        E.hydrogen : (0.85, 1.09085413633, 1.425952),
-        E.carbon : (0.72, 0.48435382330, 1.058554),
-        E.nitrogen : (0.79, 0.700147318409, 0.733599),
-        E.oxygen : (0.85, 1.06557401132, 1.061039),
-        E.fluorine : (0.88, 0.5, 0.5),
-        E.phosphorus : (0.86, 0.5, 0.5),
-        E.sulfur : (0.96, 0.602256336067, -0.703469),
-        None : (0.8, 0.5, 0.5) # default
+_SCREEN_PARAMETERS = { # normal, GBn, GBn2, GBn2 nucleic
+        E.hydrogen : (0.85, 1.09085413633, 1.425952, 1.696538),
+        E.carbon : (0.72, 0.48435382330, 1.058554, 1.268902),
+        E.nitrogen : (0.79, 0.700147318409, 0.733599, 1.4259728),
+        E.oxygen : (0.85, 1.06557401132, 1.061039, 0.1840098),
+        E.fluorine : (0.88, 0.5, 0.5, 0.5),
+        E.phosphorus : (0.86, 0.5, 0.5, 1.5450597),
+        E.sulfur : (0.96, 0.602256336067, -0.703469, 0.05),
+        None : (0.8, 0.5, 0.5, 0.5) # default
 }
 _SCREEN_PARAMETERS[E.deuterium] = _SCREEN_PARAMETERS[E.hydrogen]
 
+_NUCLEIC_ACID_RESIDUES = ['A', 'C', 'G', 'U', 'DA', 'DC', 'DG', 'DT']
+
 def _screen_parameter(atom):
     return _SCREEN_PARAMETERS.get(atom.element, _SCREEN_PARAMETERS[None])
 
@@ -873,7 +875,16 @@ class GBSAGBn2Force(GBSAGBnForce):
         E.oxygen:       [0.867814, 0.876635, 0.387882],
         E.sulfur:       [0.867814, 0.876635, 0.387882],
     }
-    _default_atom_params = [0.8, 4.85, 0.5]
+    _atom_params_nucleic = {
+        E.hydrogen:     [0.537050, 0.362861, 0.116704],
+        E.deuterium:    [0.537050, 0.362861, 0.116704],
+        E.carbon:       [0.331670, 0.196842, 0.093422],
+        E.nitrogen:     [0.686311, 0.463189, 0.138722],
+        E.oxygen:       [0.606344, 0.463006, 0.142262],
+        E.sulfur:       [0.606344, 0.463006, 0.142262],
+        E.phosphorus:   [0.418365, 0.290054, 0.1064245],
+    }
+    _default_atom_params = [1.0, 0.8, 4.851]
 
     @classmethod
     def getStandardParameters(cls, topology):
@@ -897,14 +908,23 @@ class GBSAGBn2Force(GBSAGBnForce):
             radii = numpy.empty([natoms,5], numpy.double)
             radii[:,0] = _mbondi3_radii(topology)/10
             for atom, rad in zip(topology.atoms(), radii):
-                rad[1] = _screen_parameter(atom)[2]
-                rad[2:] = cls._atom_params.get(atom.element, cls._default_atom_params)
+                if atom.residue.name in _NUCLEIC_ACID_RESIDUES:
+                    rad[1] = _screen_parameter(atom)[3]
+                    rad[2:] = cls._atom_params_nucleic.get(atom.element, cls._default_atom_params)
+                else:
+                    rad[1] = _screen_parameter(atom)[2]
+                    rad[2:] = cls._atom_params.get(atom.element, cls._default_atom_params)
         else:
             radii = [[r/10, 0, 0, 0, 0] for r in _mbondi3_radii(topology)]
             for atom, rad in zip(topology.atoms(), radii):
-                rad[1] = _screen_parameter(atom)[2]
-                for i, p in enumerate(cls._atom_params.get(atom.element, cls._default_atom_params)):
-                    rad[2+i] = p
+                if atom.residue.name in _NUCLEIC_ACID_RESIDUES:
+                    rad[1] = _screen_parameter(atom)[3]
+                    for i, p in enumerate(cls._atom_params_nucleic.get(atom.element, cls._default_atom_params)):
+                        rad[2+i] = p
+                else:
+                    rad[1] = _screen_parameter(atom)[2]
+                    for i, p in enumerate(cls._atom_params.get(atom.element, cls._default_atom_params)):
+                        rad[2+i] = p
         return radii
 
     def _addEnergyTerms(self):

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -384,5 +384,21 @@ class TestAmberPrmtopFile(unittest.TestCase):
         self.assertRaises(ValueError, lambda: f.addParticle([0, 0.9, 0.5]))
         self.assertRaises(ValueError, lambda: f.addParticle([0, 0.21, 0.5]))
 
+    def testNucleicGBParametes(self):
+        """Test that correct GB parameters are used for nucleic acids."""
+        prmtop = AmberPrmtopFile('systems/DNA_mbondi3.prmtop')
+        inpcrd = AmberInpcrdFile('systems/DNA_mbondi3.inpcrd')
+        sanderEnergy = [-19223.87993545, -19527.40433175, -19788.1070698]
+        for solvent, expectedEnergy in zip([OBC2, GBn, GBn2], sanderEnergy):
+            system = prmtop.createSystem(implicitSolvent=solvent, gbsaModel=None)
+            for f in system.getForces():
+                if isinstance(f, CustomGBForce) or isinstance(f, GBSAOBCForce):
+                    f.setForceGroup(1)
+            integrator = VerletIntegrator(0.001)
+            context = Context(system, integrator, Platform.getPlatformByName('Reference'))
+            context.setPositions(inpcrd.positions)
+            energy = context.getState(getEnergy=True, groups={1}).getPotentialEnergy().value_in_unit(kilojoules_per_mole)
+            self.assertAlmostEqual(energy, expectedEnergy, delta=5e-4*abs(energy))
+
 if __name__ == '__main__':
     unittest.main()