Skip to content

GitLab

Commit 94944787 authored by peastman's avatar peastman
Browse files

CustomGBForces apply the same energy offset as GBSAOBCForce

parent b245e772
Show whitespace changes
Inline Side-by-side
Showing with 66 additions and 77 deletions
wrappers/python/simtk/openmm/app/internal/amber_file_parser.py
View file @ 94944787
...@@ -793,19 +793,24 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -793,19 +793,24 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
if verbose: print "Adding GB parameters..." if verbose: print "Adding GB parameters..."
charges = prmtop.getCharges() charges = prmtop.getCharges()
symbls = None symbls = None
cutoff = None
if nonbondedMethod != 'NoCutoff':
cutoff = nonbondedCutoff
if units.is_quantity(cutoff):
cutoff = cutoff.value_in_unit(units.nanometers)
if gbmodel == 'GBn': if gbmodel == 'GBn':
symbls = prmtop.getAtomTypes() symbls = prmtop.getAtomTypes()
gb_parms = prmtop.getGBParms(symbls) gb_parms = prmtop.getGBParms(symbls)
if gbmodel == 'HCT': if gbmodel == 'HCT':
gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE') gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE', cutoff)
elif gbmodel == 'OBC1': elif gbmodel == 'OBC1':
gb = customgb.GBSAOBC1Force(solventDielectric, soluteDielectric, 'ACE') gb = customgb.GBSAOBC1Force(solventDielectric, soluteDielectric, 'ACE', cutoff)
elif gbmodel == 'OBC2': elif gbmodel == 'OBC2':
gb = mm.GBSAOBCForce() gb = mm.GBSAOBCForce()
gb.setSoluteDielectric(soluteDielectric) gb.setSoluteDielectric(soluteDielectric)
gb.setSolventDielectric(solventDielectric) gb.setSolventDielectric(solventDielectric)
elif gbmodel == 'GBn': elif gbmodel == 'GBn':
gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, 'ACE') gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, 'ACE', cutoff)
else: else:
raise Exception("Illegal value specified for implicit solvent model") raise Exception("Illegal value specified for implicit solvent model")
for iAtom in range(prmtop.getNumAtoms()): for iAtom in range(prmtop.getNumAtoms()):
......
wrappers/python/simtk/openmm/app/internal/customgbforces.py
View file @ 94944787
...@@ -188,12 +188,27 @@ for i in range (len(d0)): ...@@ -188,12 +188,27 @@ for i in range (len(d0)):
m0[i]=m0[i]*10 m0[i]=m0[i]*10
def _createEnergyTerms(force, SA, cutoff):
# Add the energy terms to the CustomGBForce. These are identical for all the GB models.
force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
if SA=='ACE':
force.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
elif SA is not None:
raise ValueError('Unknown surface area method: '+SA)
if cutoff is None:
force.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
else:
force.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2*(1/"+str(cutoff)+"-1/f);"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
""" """
Amber Equivalent: igb = 1 Amber Equivalent: igb = 1
""" """
def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None): def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None):
custom = CustomGBForce() custom = CustomGBForce()
...@@ -212,22 +227,13 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None): ...@@ -212,22 +227,13 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None):
custom.addComputedValue("B", "1/(1/or-I);" custom.addComputedValue("B", "1/(1/or-I);"
"or=radius-offset", CustomGBForce.SingleParticle) "or=radius-offset", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, SA, cutoff)
custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
if SA=='ACE':
custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
elif SA is not None:
raise ValueError('Unknown surface area method: '+SA)
custom.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
return custom return custom
""" """
Amber Equivalents: igb = 2 Amber Equivalents: igb = 2
""" """
def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None): def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None):
custom = CustomGBForce() custom = CustomGBForce()
...@@ -246,21 +252,13 @@ def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None): ...@@ -246,21 +252,13 @@ def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None):
custom.addComputedValue("B", "1/(1/or-tanh(0.8*psi+2.909125*psi^3)/radius);" custom.addComputedValue("B", "1/(1/or-tanh(0.8*psi+2.909125*psi^3)/radius);"
"psi=I*or; or=radius-offset", CustomGBForce.SingleParticle) "psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, SA, cutoff)
custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
if SA=='ACE':
custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
elif SA is not None:
raise ValueError('Unknown surface area method: '+SA)
custom.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
return custom return custom
""" """
Amber Equivalents: igb = 5 Amber Equivalents: igb = 5
""" """
def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None): def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None):
custom = CustomGBForce() custom = CustomGBForce()
...@@ -279,21 +277,13 @@ def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None): ...@@ -279,21 +277,13 @@ def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None):
custom.addComputedValue("B", "1/(1/or-tanh(psi-0.8*psi^2+4.85*psi^3)/radius);" custom.addComputedValue("B", "1/(1/or-tanh(psi-0.8*psi^2+4.85*psi^3)/radius);"
"psi=I*or; or=radius-offset", CustomGBForce.SingleParticle) "psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, SA, cutoff)
custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
if SA=='ACE':
custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
elif SA is not None:
raise ValueError('Unknown surface area method: '+SA)
custom.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
return custom return custom
""" """
Amber Equivalents: igb = 7 Amber Equivalents: igb = 7
""" """
def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None): def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None):
""" """
...@@ -331,13 +321,5 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None): ...@@ -331,13 +321,5 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None):
custom.addComputedValue("B", "1/(1/or-tanh(1.09511284*psi-1.907992938*psi^2+2.50798245*psi^3)/radius);" custom.addComputedValue("B", "1/(1/or-tanh(1.09511284*psi-1.907992938*psi^2+2.50798245*psi^3)/radius);"
"psi=I*or; or=radius-offset", CustomGBForce.SingleParticle) "psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, SA, cutoff)
custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
if SA=='ACE':
custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
elif SA is not None:
raise ValueError('Unknown surface area method: '+SA)
custom.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
return custom return custom
wrappers/python/tests/TestAmberPrmtopFile.py
View file @ 94944787
...@@ -93,19 +93,19 @@ class TestAmberPrmtopFile(unittest.TestCase): ...@@ -93,19 +93,19 @@ class TestAmberPrmtopFile(unittest.TestCase):
self.assertTrue(any(isinstance(f, force_type) for f in forces)) self.assertTrue(any(isinstance(f, force_type) for f in forces))
def test_ImplicitSolventParameters(self): def test_ImplicitSolventParameters(self):
"""Test that solventDielectric and soluteDielectric are passed correctly """Test that parameters are set correctly for the different types of implicit solvent."""
for the different types of implicit solvent. methodMap = {NoCutoff:NonbondedForce.NoCutoff,
CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic}
"""
for implicitSolvent_value in [HCT, OBC1, OBC2, GBn]: for implicitSolvent_value in [HCT, OBC1, OBC2, GBn]:
for method in methodMap:
system = self.prmtop2.createSystem(implicitSolvent=implicitSolvent_value, system = self.prmtop2.createSystem(implicitSolvent=implicitSolvent_value,
solventDielectric=50.0, solventDielectric=50.0, soluteDielectric=0.9, nonbondedMethod=method)
soluteDielectric = 0.9)
found_matching_solvent_dielectric=False found_matching_solvent_dielectric=False
found_matching_solute_dielectric=False found_matching_solute_dielectric=False
if implicitSolvent_value in set([HCT, OBC1, GBn]): if implicitSolvent_value in set([HCT, OBC1, GBn]):
for force in system.getForces(): for force in system.getForces():
if isinstance(force, CustomGBForce): if isinstance(force, CustomGBForce):
self.assertEqual(force.getNonbondedMethod(), methodMap[method])
for j in range(force.getNumGlobalParameters()): for j in range(force.getNumGlobalParameters()):
if (force.getGlobalParameterName(j) == 'solventDielectric' and if (force.getGlobalParameterName(j) == 'solventDielectric' and
force.getGlobalParameterDefaultValue(j) == 50.0): force.getGlobalParameterDefaultValue(j) == 50.0):
...@@ -115,20 +115,22 @@ class TestAmberPrmtopFile(unittest.TestCase): ...@@ -115,20 +115,22 @@ class TestAmberPrmtopFile(unittest.TestCase):
found_matching_solute_dielectric = True found_matching_solute_dielectric = True
if isinstance(force, NonbondedForce): if isinstance(force, NonbondedForce):
self.assertEqual(force.getReactionFieldDielectric(), 1.0) self.assertEqual(force.getReactionFieldDielectric(), 1.0)
self.assertEqual(force.getNonbondedMethod(), methodMap[method])
self.assertTrue(found_matching_solvent_dielectric and self.assertTrue(found_matching_solvent_dielectric and
found_matching_solute_dielectric) found_matching_solute_dielectric)
else: else:
for force in system.getForces(): for force in system.getForces():
if isinstance(force, GBSAOBCForce): if isinstance(force, GBSAOBCForce):
self.assertEqual(force.getNonbondedMethod(), methodMap[method])
if force.getSolventDielectric() == 50.0: if force.getSolventDielectric() == 50.0:
found_matching_solvent_dielectric = True found_matching_solvent_dielectric = True
if force.getSoluteDielectric() == 0.9: if force.getSoluteDielectric() == 0.9:
found_matching_solute_dielectric = True found_matching_solute_dielectric = True
if isinstance(force, NonbondedForce): if isinstance(force, NonbondedForce):
self.assertEqual(force.getReactionFieldDielectric(), 1.0) self.assertEqual(force.getReactionFieldDielectric(), 1.0)
self.assertEqual(force.getNonbondedMethod(), methodMap[method])
self.assertTrue(found_matching_solvent_dielectric and self.assertTrue(found_matching_solvent_dielectric and
found_matching_solute_dielectric) found_matching_solute_dielectric)
if __name__ == '__main__': if __name__ == '__main__':
unittest.main() unittest.main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment