Skip to content

GitLab

Commit 3e0925fe authored by Peter Eastman's avatar Peter Eastman
Browse files

Bug fixes

parent 8316aba0
Show whitespace changes
Inline Side-by-side
Showing with 25 additions and 25 deletions
openmmapi/include/openmm/CustomGBForce.h
View file @ 3e0925fe
...@@ -88,13 +88,13 @@ namespace OpenMM { ...@@ -88,13 +88,13 @@ namespace OpenMM {
* custom->addPerParticleParameter("scale"); * custom->addPerParticleParameter("scale");
* custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric()); * custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric());
* custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric()); * custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric());
* custom->addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r*(1/U^2-1/L^2))+0.5*log(L/U)/r+0.25*(sr2*sr2/r)*(1/L^2-1/U^2)+C);" * custom->addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
* "U=r+sr2;" * "U=r+sr2;"
* "C=2*(1/or1-1/L)*step(sr2-r-or1);" * "C=2*(1/or1-1/L)*step(sr2-r-or1);"
* "L=step(or1-D)*or1+step(D-or1)*D;" * "L=step(or1-D)*or1+(1-step(or1-D))*D;"
* "D=step(r-sr2)*(r-sr2)+step(sr2-r)*(sr2-r);" * "D=step(r-sr2)*(r-sr2)+(1-step(r-sr2))*(sr2-r);"
* "sr1 = scale1*or1; sr2 = scale2*or2;" * "sr2 = scale2*or2;"
* "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePair); * "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions);
* custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);" * custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
* "psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle); * "psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle);
* custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B", * custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B",
......
platforms/opencl/src/kernels/gbsaObc_default.cl
View file @ 3e0925fe
...@@ -130,9 +130,9 @@ __kernel void computeBornSum(__global float* global_bornSum, __local float* loca ...@@ -130,9 +130,9 @@ __kernel void computeBornSum(__global float* global_bornSum, __local float* loca
bornSum += 2.0f*(1.0f/params1.x-l_ij); bornSum += 2.0f*(1.0f/params1.x-l_ij);
} }
float rScaledRadiusI = r+params1.y; float rScaledRadiusI = r+params1.y;
if (params2.x < rScaledRadiusJ) { if (params2.x < rScaledRadiusI) {
float l_ij = 1.0f/max(params2.x, fabs(r-params1.y)); float l_ij = 1.0f/max(params2.x, fabs(r-params1.y));
float u_ij = 1.0f/rScaledRadiusJ; float u_ij = 1.0f/rScaledRadiusI;
float l_ij2 = l_ij*l_ij; float l_ij2 = l_ij*l_ij;
float u_ij2 = u_ij*u_ij; float u_ij2 = u_ij*u_ij;
float ratio = log(u_ij / l_ij); float ratio = log(u_ij / l_ij);
......
platforms/opencl/src/kernels/gbsaObc_nvidia.cl
View file @ 3e0925fe
...@@ -127,9 +127,9 @@ __kernel void computeBornSum(__global float* global_bornSum, __local float* loca ...@@ -127,9 +127,9 @@ __kernel void computeBornSum(__global float* global_bornSum, __local float* loca
bornSum += 2.0f*(1.0f/params1.x-l_ij); bornSum += 2.0f*(1.0f/params1.x-l_ij);
} }
float rScaledRadiusI = r+params1.y; float rScaledRadiusI = r+params1.y;
if (params2.x < rScaledRadiusJ) { if (params2.x < rScaledRadiusI) {
float l_ij = 1.0f/max(params2.x, fabs(r-params1.y)); float l_ij = 1.0f/max(params2.x, fabs(r-params1.y));
float u_ij = 1.0f/rScaledRadiusJ; float u_ij = 1.0f/rScaledRadiusI;
float l_ij2 = l_ij*l_ij; float l_ij2 = l_ij*l_ij;
float u_ij2 = u_ij*u_ij; float u_ij2 = u_ij*u_ij;
float ratio = log(u_ij / l_ij); float ratio = log(u_ij / l_ij);
...@@ -195,9 +195,9 @@ __kernel void computeBornSum(__global float* global_bornSum, __local float* loca ...@@ -195,9 +195,9 @@ __kernel void computeBornSum(__global float* global_bornSum, __local float* loca
bornSum += 2.0f*(1.0f/params1.x-l_ij); bornSum += 2.0f*(1.0f/params1.x-l_ij);
} }
float rScaledRadiusI = r+params1.y; float rScaledRadiusI = r+params1.y;
if (params2.x < rScaledRadiusJ) { if (params2.x < rScaledRadiusI) {
float l_ij = 1.0f/max(params2.x, fabs(r-params1.y)); float l_ij = 1.0f/max(params2.x, fabs(r-params1.y));
float u_ij = 1.0f/rScaledRadiusJ; float u_ij = 1.0f/rScaledRadiusI;
float l_ij2 = l_ij*l_ij; float l_ij2 = l_ij*l_ij;
float u_ij2 = u_ij*u_ij; float u_ij2 = u_ij*u_ij;
float ratio = log(u_ij / l_ij); float ratio = log(u_ij / l_ij);
......
platforms/opencl/tests/TestOpenCLCustomGBForce.cpp
View file @ 3e0925fe
...@@ -75,17 +75,17 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe ...@@ -75,17 +75,17 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
custom->addPerParticleParameter("scale"); custom->addPerParticleParameter("scale");
custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric()); custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric());
custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric()); custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric());
custom->addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r*(1/U^2-1/L^2))+0.5*log(L/U)/r+0.25*(sr2*sr2/r)*(1/L^2-1/U^2)+C);" custom->addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
"U=r+sr2;" "U=r+sr2;"
"C=2*(1/or1-1/L)*step(sr2-r-or1);" "C=2*(1/or1-1/L)*step(sr2-r-or1);"
"L=step(or1-D)*or1+step(D-or1)*D;" "L=step(or1-D)*or1+(1-step(or1-D))*D;"
"D=step(r-sr2)*(r-sr2)+step(sr2-r)*(sr2-r);" "D=step(r-sr2)*(r-sr2)+(1-step(r-sr2))*(sr2-r);"
"sr1 = scale1*or1; sr2 = scale2*or2;" "sr2 = scale2*or2;"
"or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions); "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions);
custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);" custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
"psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle); "psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle);
custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce::SingleParticle); custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce::SingleParticle);
custom->addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;" custom->addEnergyTerm("-138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions); "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions);
vector<Vec3> positions(numParticles); vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles); vector<Vec3> velocities(numParticles);
...@@ -175,10 +175,10 @@ void testTabulatedFunction(bool interpolating) { ...@@ -175,10 +175,10 @@ void testTabulatedFunction(bool interpolating) {
int main() { int main() {
try { try {
testOBC(GBSAOBCForce::NoCutoff, CustomGBForce::NoCutoff); testOBC(GBSAOBCForce::NoCutoff, CustomGBForce::NoCutoff);
// testOBC(GBSAOBCForce::CutoffNonPeriodic, CustomGBForce::CutoffNonPeriodic); testOBC(GBSAOBCForce::CutoffNonPeriodic, CustomGBForce::CutoffNonPeriodic);
// testOBC(GBSAOBCForce::CutoffPeriodic, CustomGBForce::CutoffPeriodic); testOBC(GBSAOBCForce::CutoffPeriodic, CustomGBForce::CutoffPeriodic);
// testTabulatedFunction(true); testTabulatedFunction(true);
// testTabulatedFunction(false); testTabulatedFunction(false);
} }
catch(const exception& e) { catch(const exception& e) {
cout << "exception: " << e.what() << endl; cout << "exception: " << e.what() << endl;
......
platforms/reference/tests/TestReferenceCustomGBForce.cpp
View file @ 3e0925fe
...@@ -75,12 +75,12 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe ...@@ -75,12 +75,12 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
custom->addPerParticleParameter("scale"); custom->addPerParticleParameter("scale");
custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric()); custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric());
custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric()); custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric());
custom->addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r*(1/U^2-1/L^2))+0.5*log(L/U)/r+0.25*(sr2*sr2/r)*(1/L^2-1/U^2)+C);" custom->addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
"U=r+sr2;" "U=r+sr2;"
"C=2*(1/or1-1/L)*step(sr2-r-or1);" "C=2*(1/or1-1/L)*step(sr2-r-or1);"
"L=step(or1-D)*or1+step(D-or1)*D;" "L=step(or1-D)*or1+(1-step(or1-D))*D;"
"D=step(r-sr2)*(r-sr2)+step(sr2-r)*(sr2-r);" "D=step(r-sr2)*(r-sr2)+(1-step(r-sr2))*(sr2-r);"
"sr1 = scale1*or1; sr2 = scale2*or2;" "sr2 = scale2*or2;"
"or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions); "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions);
custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);" custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
"psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle); "psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle);
......
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