Skip to content

GitLab

Commit 71b67421 authored by Peter Eastman's avatar Peter Eastman
Browse files

Further optimizations to custom GB models

parent 92a4bbe1
Hide whitespace changes
Inline Side-by-side
Showing with 16 additions and 21 deletions
wrappers/python/simtk/openmm/app/internal/customgbforces.py
View file @ 71b67421
...@@ -189,10 +189,12 @@ for i in range (len(d0)): ...@@ -189,10 +189,12 @@ for i in range (len(d0)):
m0[i]=m0[i]*10 m0[i]=m0[i]*10
def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, kappa): def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, kappa, offset):
# Add the energy terms to the CustomGBForce. These are identical for all the GB models. # Add the energy terms to the CustomGBForce. These are identical for all the GB models.
params = "; solventDielectric="+str(solventDielectric)+"; soluteDielectric="+str(soluteDielectric)+"; kappa="+str(kappa)+"; cutoff="+str(cutoff) params = "; solventDielectric=%.16g; soluteDielectric=%.16g; kappa=%.16g; offset=%.16g" % (solventDielectric, soluteDielectric, kappa, offset)
if cutoff is not None:
params += "; cutoff=%.16g" % cutoff
if kappa > 0: if kappa > 0:
force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-exp(-kappa*B)/solventDielectric)*q^2/B"+params, force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-exp(-kappa*B)/solventDielectric)*q^2/B"+params,
CustomGBForce.SingleParticle) CustomGBForce.SingleParticle)
...@@ -203,7 +205,7 @@ def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, k ...@@ -203,7 +205,7 @@ def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, k
force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B"+params, force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B"+params,
CustomGBForce.SingleParticle) CustomGBForce.SingleParticle)
if SA=='ACE': if SA=='ACE':
force.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6; radius=or+offset", CustomGBForce.SingleParticle) force.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6; radius=or+offset"+params, CustomGBForce.SingleParticle)
elif SA is not None: elif SA is not None:
raise ValueError('Unknown surface area method: '+SA) raise ValueError('Unknown surface area method: '+SA)
if cutoff is None: if cutoff is None:
...@@ -234,14 +236,13 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None, ...@@ -234,14 +236,13 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom.addPerParticleParameter("q") custom.addPerParticleParameter("q")
custom.addPerParticleParameter("or") # Offset radius custom.addPerParticleParameter("or") # Offset radius
custom.addPerParticleParameter("sr") # Scaled offset radius custom.addPerParticleParameter("sr") # Scaled offset radius
custom.addGlobalParameter("offset", 0.009)
custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);" custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
"U=r+sr2;" "U=r+sr2;"
"L=max(or1, D);" "L=max(or1, D);"
"D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions) "D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
custom.addComputedValue("B", "1/(1/or-I)", CustomGBForce.SingleParticle) custom.addComputedValue("B", "1/(1/or-I)", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa) _createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
return custom return custom
""" """
...@@ -255,15 +256,14 @@ def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None, ...@@ -255,15 +256,14 @@ def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom.addPerParticleParameter("q") custom.addPerParticleParameter("q")
custom.addPerParticleParameter("or") # Offset radius custom.addPerParticleParameter("or") # Offset radius
custom.addPerParticleParameter("sr") # Scaled offset radius custom.addPerParticleParameter("sr") # Scaled offset radius
custom.addGlobalParameter("offset", 0.009)
custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);" custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
"U=r+sr2;" "U=r+sr2;"
"L=max(or1, D);" "L=max(or1, D);"
"D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions) "D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
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; radius=or+offset", CustomGBForce.SingleParticle) "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa) _createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
return custom return custom
""" """
...@@ -277,15 +277,14 @@ def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, ...@@ -277,15 +277,14 @@ def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom.addPerParticleParameter("q") custom.addPerParticleParameter("q")
custom.addPerParticleParameter("or") # Offset radius custom.addPerParticleParameter("or") # Offset radius
custom.addPerParticleParameter("sr") # Scaled offset radius custom.addPerParticleParameter("sr") # Scaled offset radius
custom.addGlobalParameter("offset", 0.009)
custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);" custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
"U=r+sr2;" "U=r+sr2;"
"L=max(or1, D);" "L=max(or1, D);"
"D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions) "D=abs(r-sr2)", CustomGBForce.ParticlePairNoExclusions)
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; radius=or+offset", CustomGBForce.SingleParticle) "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa) _createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
return custom return custom
""" """
...@@ -309,8 +308,6 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None, ...@@ -309,8 +308,6 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom.addPerParticleParameter("or") # Offset radius custom.addPerParticleParameter("or") # Offset radius
custom.addPerParticleParameter("sr") # Scaled offset radius custom.addPerParticleParameter("sr") # Scaled offset radius
custom.addGlobalParameter("offset", 0.009)
custom.addTabulatedFunction("getd0", Discrete1DFunction(d0)) custom.addTabulatedFunction("getd0", Discrete1DFunction(d0))
custom.addTabulatedFunction("getm0", Discrete1DFunction(m0)) custom.addTabulatedFunction("getm0", Discrete1DFunction(m0))
...@@ -322,11 +319,11 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None, ...@@ -322,11 +319,11 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None,
"L=max(or1, D);" "L=max(or1, D);"
"D=abs(r-sr2);" "D=abs(r-sr2);"
"radius1=or1+offset; radius2=or2+offset;" "radius1=or1+offset; radius2=or2+offset;"
"neckScale=0.361825; neckCut=0.68", CustomGBForce.ParticlePairNoExclusions) "neckScale=0.361825; neckCut=0.68; offset=0.009", CustomGBForce.ParticlePairNoExclusions)
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; radius=or+offset", CustomGBForce.SingleParticle) "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa) _createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
return custom return custom
""" """
...@@ -353,8 +350,6 @@ def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, ...@@ -353,8 +350,6 @@ def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom.addPerParticleParameter("beta") custom.addPerParticleParameter("beta")
custom.addPerParticleParameter("gamma") custom.addPerParticleParameter("gamma")
custom.addGlobalParameter("offset", 0.0195141)
custom.addTabulatedFunction("getd0", Discrete1DFunction(d0)) custom.addTabulatedFunction("getd0", Discrete1DFunction(d0))
custom.addTabulatedFunction("getm0", Discrete1DFunction(m0)) custom.addTabulatedFunction("getm0", Discrete1DFunction(m0))
...@@ -366,11 +361,11 @@ def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, ...@@ -366,11 +361,11 @@ def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
"L=max(or1, D);" "L=max(or1, D);"
"D=abs(r-sr2);" "D=abs(r-sr2);"
"radius1=or1+offset; radius2=or2+offset;" "radius1=or1+offset; radius2=or2+offset;"
"neckScale=0.826836; neckCut=0.68", CustomGBForce.ParticlePairNoExclusions) "neckScale=0.826836; neckCut=0.68; offset=0.0195141", CustomGBForce.ParticlePairNoExclusions)
custom.addComputedValue("B", "1/(1/or-tanh(alpha*psi-beta*psi^2+gamma*psi^3)/radius);" custom.addComputedValue("B", "1/(1/or-tanh(alpha*psi-beta*psi^2+gamma*psi^3)/radius);"
"psi=I*or; radius=or+offset", CustomGBForce.SingleParticle) "psi=I*or; radius=or+offset; offset=0.0195141", CustomGBForce.SingleParticle)
_createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa) _createEnergyTerms(custom, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.0195141)
return custom return custom
......
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