Commit 4c450bf5 authored by Jason Swails's avatar Jason Swails
Browse files

Move to using a continuous 2D function.

parent dafd0f97
......@@ -31,7 +31,7 @@ USE OR OTHER DEALINGS IN THE SOFTWARE.
from __future__ import division
from simtk.openmm import CustomGBForce, Discrete2DFunction
from simtk.openmm import CustomGBForce, Continuous2DFunction
d0=[2.26685,2.32548,2.38397,2.44235,2.50057,2.55867,2.61663,2.67444,
2.73212,2.78965,2.84705,2.9043,2.96141,3.0184,3.07524,3.13196,
......@@ -308,12 +308,11 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom.addPerParticleParameter("or") # Offset radius
custom.addPerParticleParameter("sr") # Scaled offset radius
custom.addTabulatedFunction("getd0", Discrete2DFunction(21, 21, d0))
custom.addTabulatedFunction("getm0", Discrete2DFunction(21, 21, m0))
custom.addTabulatedFunction("getd0", Continuous2DFunction(21, 21, d0, 0.1, 0.2, 0.1, 0.2))
custom.addTabulatedFunction("getm0", Continuous2DFunction(21, 21, m0, 0.1, 0.2, 0.1, 0.2))
custom.addComputedValue("I", "Ivdw+neckScale*Ineck;"
"Ineck=step(radius1+radius2+neckCut-r)*getm0(index1,index2)/(1+100*(r-getd0(index1,index2))^2+0.3*1000000*(r-getd0(index1,index2))^6);"
"index1=radius1*200-20; index2=radius2*200-20;"
"Ineck=step(radius1+radius2+neckCut-r)*getm0(radius1,radius2)/(1+100*(r-getd0(radius1,radius2))^2+0.3*1000000*(r-getd0(radius1,radius2))^6);"
"Ivdw=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;"
"L=max(or1, D);"
......@@ -350,12 +349,11 @@ def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
custom.addPerParticleParameter("beta")
custom.addPerParticleParameter("gamma")
custom.addTabulatedFunction("getd0", Discrete2DFunction(21, 21, d0))
custom.addTabulatedFunction("getm0", Discrete2DFunction(21, 21, m0))
custom.addTabulatedFunction("getd0", Continuous2DFunction(21, 21, d0, 0.1, 0.2, 0.1, 0.2))
custom.addTabulatedFunction("getm0", Continuous2DFunction(21, 21, m0, 0.1, 0.2, 0.1, 0.2))
custom.addComputedValue("I", "Ivdw+neckScale*Ineck;"
"Ineck=step(radius1+radius2+neckCut-r)*getm0(index1,index2)/(1+100*(r-getd0(index1,index2))^2+0.3*1000000*(r-getd0(index1,index2))^6);"
"index1=radius1*200-20; index2=radius2*200-20;"
"Ineck=step(radius1+radius2+neckCut-r)*getm0(radius1,radius2)/(1+100*(r-getd0(radius1,radius2))^2+0.3*1000000*(r-getd0(radius1,radius2))^6);"
"Ivdw=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;"
"L=max(or1, D);"
......
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