Skip to content

GitLab

Commit 9df60d18 authored by peastman's avatar peastman
Browse files

Merge pull request #123 from peastman/gb2 

Apply an offset to GB energy when using a cutoff to reduce the discontinuity
parents 1d4f5411 f55a30ae
Show whitespace changes
Inline Side-by-side
Showing with 165 additions and 121 deletions
platforms/cuda/src/kernels/gbsaObc1.cu
View file @ 9df60d18
......@@ -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 @ 9df60d18
......@@ -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 @ 9df60d18
......@@ -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 @ 9df60d18
......@@ -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 @ 9df60d18
......@@ -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 @ 9df60d18
......@@ -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 @ 9df60d18
......@@ -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 @ 9df60d18
......@@ -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);
......
wrappers/python/simtk/openmm/app/internal/amber_file_parser.py
View file @ 9df60d18
......@@ -793,19 +793,24 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
if verbose: print "Adding GB parameters..."
charges = prmtop.getCharges()
symbls = None
cutoff = None
if nonbondedMethod != 'NoCutoff':
cutoff = nonbondedCutoff
if units.is_quantity(cutoff):
cutoff = cutoff.value_in_unit(units.nanometers)
if gbmodel == 'GBn':
symbls = prmtop.getAtomTypes()
gb_parms = prmtop.getGBParms(symbls)
if gbmodel == 'HCT':
gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE')
gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE', cutoff)
elif gbmodel == 'OBC1':
gb = customgb.GBSAOBC1Force(solventDielectric, soluteDielectric, 'ACE')
gb = customgb.GBSAOBC1Force(solventDielectric, soluteDielectric, 'ACE', cutoff)
elif gbmodel == 'OBC2':
gb = mm.GBSAOBCForce()
gb.setSoluteDielectric(soluteDielectric)
gb.setSolventDielectric(solventDielectric)
elif gbmodel == 'GBn':
gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, 'ACE')
gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, 'ACE', cutoff)
else:
raise Exception("Illegal value specified for implicit solvent model")
for iAtom in range(prmtop.getNumAtoms()):
......
wrappers/python/simtk/openmm/app/internal/customgbforces.py
View file @ 9df60d18
......@@ -188,12 +188,27 @@ for i in range (len(d0)):
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/f-"+str(1/cutoff)+");"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
"""
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()
......@@ -212,22 +227,13 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None):
custom.addComputedValue("B", "1/(1/or-I);"
"or=radius-offset", CustomGBForce.SingleParticle)
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)
_createEnergyTerms(custom, SA, cutoff)
return custom
"""
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()
......@@ -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);"
"psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
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)
_createEnergyTerms(custom, SA, cutoff)
return custom
"""
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()
......@@ -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);"
"psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
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)
_createEnergyTerms(custom, SA, cutoff)
return custom
"""
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):
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)
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)
_createEnergyTerms(custom, SA, cutoff)
return custom
wrappers/python/tests/TestAmberPrmtopFile.py
View file @ 9df60d18
......@@ -93,19 +93,19 @@ class TestAmberPrmtopFile(unittest.TestCase):
self.assertTrue(any(isinstance(f, force_type) for f in forces))
def test_ImplicitSolventParameters(self):
"""Test that solventDielectric and soluteDielectric are passed correctly
for the different types of implicit solvent.
"""
"""Test that parameters are set correctly for the different types of implicit solvent."""
methodMap = {NoCutoff:NonbondedForce.NoCutoff,
CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic}
for implicitSolvent_value in [HCT, OBC1, OBC2, GBn]:
for method in methodMap:
system = self.prmtop2.createSystem(implicitSolvent=implicitSolvent_value,
solventDielectric=50.0,
soluteDielectric = 0.9)
solventDielectric=50.0, soluteDielectric=0.9, nonbondedMethod=method)
found_matching_solvent_dielectric=False
found_matching_solute_dielectric=False
if implicitSolvent_value in set([HCT, OBC1, GBn]):
for force in system.getForces():
if isinstance(force, CustomGBForce):
self.assertEqual(force.getNonbondedMethod(), methodMap[method])
for j in range(force.getNumGlobalParameters()):
if (force.getGlobalParameterName(j) == 'solventDielectric' and
force.getGlobalParameterDefaultValue(j) == 50.0):
......@@ -115,20 +115,22 @@ class TestAmberPrmtopFile(unittest.TestCase):
found_matching_solute_dielectric = True
if isinstance(force, NonbondedForce):
self.assertEqual(force.getReactionFieldDielectric(), 1.0)
self.assertEqual(force.getNonbondedMethod(), methodMap[method])
self.assertTrue(found_matching_solvent_dielectric and
found_matching_solute_dielectric)
else:
for force in system.getForces():
if isinstance(force, GBSAOBCForce):
self.assertEqual(force.getNonbondedMethod(), methodMap[method])
if force.getSolventDielectric() == 50.0:
found_matching_solvent_dielectric = True
if force.getSoluteDielectric() == 0.9:
found_matching_solute_dielectric = True
if isinstance(force, NonbondedForce):
self.assertEqual(force.getReactionFieldDielectric(), 1.0)
self.assertEqual(force.getNonbondedMethod(), methodMap[method])
self.assertTrue(found_matching_solvent_dielectric and
found_matching_solute_dielectric)
if __name__ == '__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