Skip to content

GitLab

Commit b245e772 authored by peastman's avatar peastman
Browse files

GBSAOBCForce adds an offset to the energy to reduce the discontinuity at the cutoff

parent e59e6f93
Show whitespace changes
Inline Side-by-side
Showing with 99 additions and 44 deletions
platforms/cuda/src/kernels/gbsaObc1.cu
View file @ b245e772
......@@ -468,11 +468,16 @@ extern "C" __global__ void computeGBSAForce1(unsigned long long* __restrict__ fo
real expTerm = EXP(-D_ij);
real denominator2 = r2 + alpha2_ij*expTerm;
real denominator = SQRT(denominator2);
real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator);
real scaledChargeProduct = PREFACTOR*posq1.w*posq2.w;
real tempEnergy = scaledChargeProduct*RECIP(denominator);
real Gpol = tempEnergy*RECIP(denominator2);
real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij);
real dEdR = Gpol*(1.0f - 0.25f*expTerm);
force.w += dGpol_dalpha2_ij*bornRadius2;
#ifdef USE_CUTOFF
if (atom1 != y*TILE_SIZE+j)
tempEnergy -= scaledChargeProduct/CUTOFF;
#endif
energy += 0.5f*tempEnergy;
delta *= dEdR;
force.x -= delta.x;
......@@ -520,11 +525,15 @@ extern "C" __global__ void computeGBSAForce1(unsigned long long* __restrict__ fo
real expTerm = EXP(-D_ij);
real denominator2 = r2 + alpha2_ij*expTerm;
real denominator = SQRT(denominator2);
real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator);
real scaledChargeProduct = PREFACTOR*posq1.w*posq2.w;
real tempEnergy = scaledChargeProduct*RECIP(denominator);
real Gpol = tempEnergy*RECIP(denominator2);
real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij);
real dEdR = Gpol*(1.0f - 0.25f*expTerm);
force.w += dGpol_dalpha2_ij*bornRadius2;
#ifdef USE_CUTOFF
tempEnergy -= scaledChargeProduct/CUTOFF;
#endif
energy += tempEnergy;
delta *= dEdR;
force.x -= delta.x;
......@@ -670,11 +679,15 @@ extern "C" __global__ void computeGBSAForce1(unsigned long long* __restrict__ fo
real expTerm = EXP(-D_ij);
real denominator2 = r2 + alpha2_ij*expTerm;
real denominator = SQRT(denominator2);
real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator);
real scaledChargeProduct = PREFACTOR*posq1.w*posq2.w;
real tempEnergy = scaledChargeProduct*RECIP(denominator);
real Gpol = tempEnergy*RECIP(denominator2);
real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij);
real dEdR = Gpol*(1.0f - 0.25f*expTerm);
force.w += dGpol_dalpha2_ij*bornRadius2;
#ifdef USE_CUTOFF
tempEnergy -= scaledChargeProduct/CUTOFF;
#endif
energy += tempEnergy;
delta *= dEdR;
force.x -= delta.x;
......@@ -717,11 +730,15 @@ extern "C" __global__ void computeGBSAForce1(unsigned long long* __restrict__ fo
real expTerm = EXP(-D_ij);
real denominator2 = r2 + alpha2_ij*expTerm;
real denominator = SQRT(denominator2);
real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator);
real scaledChargeProduct = PREFACTOR*posq1.w*posq2.w;
real tempEnergy = scaledChargeProduct*RECIP(denominator);
real Gpol = tempEnergy*RECIP(denominator2);
real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij);
real dEdR = Gpol*(1.0f - 0.25f*expTerm);
force.w += dGpol_dalpha2_ij*bornRadius2;
#ifdef USE_CUTOFF
tempEnergy -= scaledChargeProduct/CUTOFF;
#endif
energy += tempEnergy;
delta *= dEdR;
force.x -= delta.x;
......
platforms/cuda/tests/TestCudaCustomGBForce.cpp
View file @ b245e772
......@@ -56,6 +56,7 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
const int numMolecules = 70;
const int numParticles = numMolecules*2;
const double boxSize = 10.0;
const double cutoff = 2.0;
// Create two systems: one with a GBSAOBCForce, and one using a CustomGBForce to implement the same interaction.
......@@ -69,8 +70,8 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
customSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
GBSAOBCForce* obc = new GBSAOBCForce();
CustomGBForce* custom = new CustomGBForce();
obc->setCutoffDistance(2.0);
custom->setCutoffDistance(2.0);
obc->setCutoffDistance(cutoff);
custom->setCutoffDistance(cutoff);
custom->addPerParticleParameter("q");
custom->addPerParticleParameter("radius");
custom->addPerParticleParameter("scale");
......@@ -86,7 +87,13 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
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);
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.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
string invCutoffString = "";
if (obcMethod != GBSAOBCForce::NoCutoff) {
stringstream s;
s<<(1.0/cutoff);
invCutoffString = s.str();
}
custom->addEnergyTerm("138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2*("+invCutoffString+"-1/f);"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
......
platforms/opencl/src/kernels/gbsaObc.cl
View file @ b245e772
......@@ -444,11 +444,16 @@ __kernel void computeGBSAForce1(
real expTerm = EXP(-D_ij);
real denominator2 = r2 + alpha2_ij*expTerm;
real denominator = SQRT(denominator2);
real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator);
real scaledChargeProduct = PREFACTOR*posq1.w*posq2.w;
real tempEnergy = scaledChargeProduct*RECIP(denominator);
real Gpol = tempEnergy*RECIP(denominator2);
real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij);
real dEdR = Gpol*(1.0f - 0.25f*expTerm);
force.w += dGpol_dalpha2_ij*bornRadius2;
#ifdef USE_CUTOFF
if (atom1 != y*TILE_SIZE+j)
tempEnergy -= scaledChargeProduct/CUTOFF;
#endif
energy += 0.5f*tempEnergy;
delta.xyz *= dEdR;
force.xyz -= delta.xyz;
......@@ -494,11 +499,15 @@ __kernel void computeGBSAForce1(
real expTerm = EXP(-D_ij);
real denominator2 = r2 + alpha2_ij*expTerm;
real denominator = SQRT(denominator2);
real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator);
real scaledChargeProduct = PREFACTOR*posq1.w*posq2.w;
real tempEnergy = scaledChargeProduct*RECIP(denominator);
real Gpol = tempEnergy*RECIP(denominator2);
real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij);
real dEdR = Gpol*(1.0f - 0.25f*expTerm);
force.w += dGpol_dalpha2_ij*bornRadius2;
#ifdef USE_CUTOFF
tempEnergy -= scaledChargeProduct/CUTOFF;
#endif
energy += tempEnergy;
delta.xyz *= dEdR;
force.xyz -= delta.xyz;
......@@ -657,11 +666,15 @@ __kernel void computeGBSAForce1(
real expTerm = EXP(-D_ij);
real denominator2 = r2 + alpha2_ij*expTerm;
real denominator = SQRT(denominator2);
real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator);
real scaledChargeProduct = PREFACTOR*posq1.w*posq2.w;
real tempEnergy = scaledChargeProduct*RECIP(denominator);
real Gpol = tempEnergy*RECIP(denominator2);
real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij);
real dEdR = Gpol*(1.0f - 0.25f*expTerm);
force.w += dGpol_dalpha2_ij*bornRadius2;
#ifdef USE_CUTOFF
tempEnergy -= scaledChargeProduct/CUTOFF;
#endif
energy += tempEnergy;
delta.xyz *= dEdR;
force.xyz -= delta.xyz;
......@@ -701,11 +714,15 @@ __kernel void computeGBSAForce1(
real expTerm = EXP(-D_ij);
real denominator2 = r2 + alpha2_ij*expTerm;
real denominator = SQRT(denominator2);
real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator);
real scaledChargeProduct = PREFACTOR*posq1.w*posq2.w;
real tempEnergy = scaledChargeProduct*RECIP(denominator);
real Gpol = tempEnergy*RECIP(denominator2);
real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij);
real dEdR = Gpol*(1.0f - 0.25f*expTerm);
force.w += dGpol_dalpha2_ij*bornRadius2;
#ifdef USE_CUTOFF
tempEnergy -= scaledChargeProduct/CUTOFF;
#endif
energy += tempEnergy;
delta.xyz *= dEdR;
force.xyz -= delta.xyz;
......
platforms/opencl/tests/TestOpenCLCustomGBForce.cpp
View file @ b245e772
......@@ -56,6 +56,7 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
const int numMolecules = 70;
const int numParticles = numMolecules*2;
const double boxSize = 10.0;
const double cutoff = 2.0;
// Create two systems: one with a GBSAOBCForce, and one using a CustomGBForce to implement the same interaction.
......@@ -69,8 +70,8 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
customSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
GBSAOBCForce* obc = new GBSAOBCForce();
CustomGBForce* custom = new CustomGBForce();
obc->setCutoffDistance(2.0);
custom->setCutoffDistance(2.0);
obc->setCutoffDistance(cutoff);
custom->setCutoffDistance(cutoff);
custom->addPerParticleParameter("q");
custom->addPerParticleParameter("radius");
custom->addPerParticleParameter("scale");
......@@ -86,7 +87,13 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
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);
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.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
string invCutoffString = "";
if (obcMethod != GBSAOBCForce::NoCutoff) {
stringstream s;
s<<(1.0/cutoff);
invCutoffString = s.str();
}
custom->addEnergyTerm("138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2*("+invCutoffString+"-1/f);"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
......
platforms/reference/include/ObcParameters.h
View file @ b245e772
......@@ -305,7 +305,7 @@ class ObcParameters {
--------------------------------------------------------------------------------------- */
bool getUseCutoff();
bool getUseCutoff() const;
/**---------------------------------------------------------------------------------------
......@@ -313,7 +313,7 @@ class ObcParameters {
--------------------------------------------------------------------------------------- */
RealOpenMM getCutoffDistance();
RealOpenMM getCutoffDistance() const;
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/gbsa/CpuObc.cpp
View file @ b245e772
......@@ -308,23 +308,18 @@ RealOpenMM CpuObc::computeBornEnergyForces( const vector<RealVec>& atomCoordinat
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
static const RealOpenMM eighth = static_cast<RealOpenMM>( 0.125 );
// ---------------------------------------------------------------------------------------
const ObcParameters* obcParameters = getObcParameters();
const int numberOfAtoms = obcParameters->getNumberOfAtoms();
// ---------------------------------------------------------------------------------------
// constants
const int numberOfAtoms = _obcParameters->getNumberOfAtoms();
const RealOpenMM dielectricOffset = _obcParameters->getDielectricOffset();
const RealOpenMM cutoffDistance = _obcParameters->getCutoffDistance();
const RealOpenMM soluteDielectric = _obcParameters->getSoluteDielectric();
const RealOpenMM solventDielectric = _obcParameters->getSolventDielectric();
RealOpenMM preFactor;
if( obcParameters->getSoluteDielectric() != zero && obcParameters->getSolventDielectric() != zero ){
preFactor = two*obcParameters->getElectricConstant()*( (one/obcParameters->getSoluteDielectric()) - (one/obcParameters->getSolventDielectric()) );
} else {
if (soluteDielectric != zero && solventDielectric != zero)
preFactor = two*_obcParameters->getElectricConstant()*((one/soluteDielectric) - (one/solventDielectric));
else
preFactor = zero;
}
const RealOpenMM dielectricOffset = obcParameters->getDielectricOffset();
// ---------------------------------------------------------------------------------------
......@@ -343,7 +338,7 @@ RealOpenMM CpuObc::computeBornEnergyForces( const vector<RealVec>& atomCoordinat
// compute the nonpolar solvation via ACE approximation
if( includeAceApproximation() ){
computeAceNonPolarForce( obcParameters, bornRadii, &obcEnergy, bornForces );
computeAceNonPolarForce( _obcParameters, bornRadii, &obcEnergy, bornForces );
}
// ---------------------------------------------------------------------------------------
......@@ -360,7 +355,7 @@ RealOpenMM CpuObc::computeBornEnergyForces( const vector<RealVec>& atomCoordinat
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcParameters->getCutoffDistance())
if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > cutoffDistance)
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
......@@ -380,8 +375,13 @@ RealOpenMM CpuObc::computeBornEnergyForces( const vector<RealVec>& atomCoordinat
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*( one + D_ij )/denominator2;
RealOpenMM energy = Gpol;
if( atomI != atomJ ){
if (_obcParameters->getUseCutoff())
energy -= partialChargeI*partialCharges[atomJ]/cutoffDistance;
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
deltaX *= dGpol_dr;
......@@ -397,10 +397,10 @@ RealOpenMM CpuObc::computeBornEnergyForces( const vector<RealVec>& atomCoordinat
inputForces[atomJ][2] -= deltaZ;
} else {
Gpol *= half;
energy *= half;
}
obcEnergy += Gpol;
obcEnergy += energy;
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
}
......@@ -411,12 +411,12 @@ RealOpenMM CpuObc::computeBornEnergyForces( const vector<RealVec>& atomCoordinat
// second main loop
const RealOpenMMVector& obcChain = getObcChain();
const RealOpenMMVector& atomicRadii = obcParameters->getAtomicRadii();
const RealOpenMMVector& atomicRadii = _obcParameters->getAtomicRadii();
const RealOpenMM alphaObc = obcParameters->getAlphaObc();
const RealOpenMM betaObc = obcParameters->getBetaObc();
const RealOpenMM gammaObc = obcParameters->getGammaObc();
const RealOpenMMVector& scaledRadiusFactor = obcParameters->getScaledRadiusFactors();
const RealOpenMM alphaObc = _obcParameters->getAlphaObc();
const RealOpenMM betaObc = _obcParameters->getBetaObc();
const RealOpenMM gammaObc = _obcParameters->getGammaObc();
const RealOpenMMVector& scaledRadiusFactor = _obcParameters->getScaledRadiusFactors();
// compute factor that depends only on the outer loop index
......@@ -440,7 +440,7 @@ RealOpenMM CpuObc::computeBornEnergyForces( const vector<RealVec>& atomCoordinat
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcParameters->getCutoffDistance())
if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > cutoffDistance)
continue;
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
......
platforms/reference/src/gbsa/ObcParameters.cpp
View file @ b245e772
......@@ -351,7 +351,7 @@ void ObcParameters::setUseCutoff( RealOpenMM distance ) {
--------------------------------------------------------------------------------------- */
bool ObcParameters::getUseCutoff() {
bool ObcParameters::getUseCutoff() const {
return _cutoff;
}
......@@ -361,7 +361,7 @@ bool ObcParameters::getUseCutoff() {
--------------------------------------------------------------------------------------- */
RealOpenMM ObcParameters::getCutoffDistance() {
RealOpenMM ObcParameters::getCutoffDistance() const {
return _cutoffDistance;
}
......
platforms/reference/tests/TestReferenceCustomGBForce.cpp
View file @ b245e772
......@@ -57,6 +57,7 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
const int numMolecules = 70;
const int numParticles = numMolecules*2;
const double boxSize = 10.0;
const double cutoff = 2.0;
ReferencePlatform platform;
// Create two systems: one with a GBSAOBCForce, and one using a CustomGBForce to implement the same interaction.
......@@ -71,8 +72,8 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
customSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
GBSAOBCForce* obc = new GBSAOBCForce();
CustomGBForce* custom = new CustomGBForce();
obc->setCutoffDistance(2.0);
custom->setCutoffDistance(2.0);
obc->setCutoffDistance(cutoff);
custom->setCutoffDistance(cutoff);
custom->addPerParticleParameter("q");
custom->addPerParticleParameter("radius");
custom->addPerParticleParameter("scale");
......@@ -88,7 +89,13 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
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);
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("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
string invCutoffString = "";
if (obcMethod != GBSAOBCForce::NoCutoff) {
stringstream s;
s<<(1.0/cutoff);
invCutoffString = s.str();
}
custom->addEnergyTerm("138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2*("+invCutoffString+"-1/f);"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
......
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