Skip to content

GitLab

Commit 36ad9ee7 authored by Michael Schnieders's avatar Michael Schnieders
Browse files

AMOEBA vdW softcore is propogated through the Reference platform; pending some... 

AMOEBA vdW softcore is propogated through the Reference platform; pending some feedback adding it to the CUDA code can begin
parent ae4920a8
Hide whitespace changes
Inline Side-by-side
Showing with 287 additions and 74 deletions
plugins/amoeba/openmmapi/include/openmm/AmoebaVdwForce.h
View file @ 36ad9ee7
......@@ -312,9 +312,9 @@ private:
NonbondedMethod nonbondedMethod;
double cutoff;
bool useDispersionCorrection;
AlchemicalMethod alchemicalMethod = AlchemicalMethod.None;
int n = 5;
double alpha = 0.7;
AlchemicalMethod alchemicalMethod;
int n;
double alpha;
std::string sigmaCombiningRule;
std::string epsilonCombiningRule;
......
plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp
View file @ 36ad9ee7
......@@ -38,28 +38,30 @@ using namespace OpenMM;
using std::string;
using std::vector;
AmoebaVdwForce::AmoebaVdwForce() : nonbondedMethod(NoCutoff), sigmaCombiningRule("CUBIC-MEAN"), epsilonCombiningRule("HHG"), cutoff(1.0e+10), useDispersionCorrection(true) {
AmoebaVdwForce::AmoebaVdwForce() : nonbondedMethod(NoCutoff), sigmaCombiningRule("CUBIC-MEAN"), epsilonCombiningRule("HHG"), cutoff(1.0e+10), useDispersionCorrection(true), alchemicalMethod(None), n(5), alpha(0.7) {
}
int AmoebaVdwForce::addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor) {
parameters.push_back(VdwInfo(parentIndex, sigma, epsilon, reductionFactor));
int AmoebaVdwForce::addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor, bool isAlchemical) {
parameters.push_back(VdwInfo(parentIndex, sigma, epsilon, reductionFactor, isAlchemical));
return parameters.size()-1;
}
void AmoebaVdwForce::getParticleParameters(int particleIndex, int& parentIndex,
double& sigma, double& epsilon, double& reductionFactor) const {
double& sigma, double& epsilon, double& reductionFactor, bool& isAlchemical) const {
parentIndex = parameters[particleIndex].parentIndex;
sigma = parameters[particleIndex].sigma;
epsilon = parameters[particleIndex].epsilon;
reductionFactor = parameters[particleIndex].reductionFactor;
isAlchemical = parameters[particleIndex].isAlchemical;
}
void AmoebaVdwForce::setParticleParameters(int particleIndex, int parentIndex,
double sigma, double epsilon, double reductionFactor) {
double sigma, double epsilon, double reductionFactor, bool isAlchemical) {
parameters[particleIndex].parentIndex = parentIndex;
parameters[particleIndex].sigma = sigma;
parameters[particleIndex].epsilon = epsilon;
parameters[particleIndex].reductionFactor = reductionFactor;
parameters[particleIndex].isAlchemical = isAlchemical;
}
void AmoebaVdwForce::setSigmaCombiningRule(const std::string& inputSigmaCombiningRule) {
......@@ -128,6 +130,33 @@ void AmoebaVdwForce::setNonbondedMethod(NonbondedMethod method) {
nonbondedMethod = method;
}
AmoebaVdwForce::AlchemicalMethod AmoebaVdwForce::getAlchemicalMethod() const {
return alchemicalMethod;
}
void AmoebaVdwForce::setAlchemicalMethod(AlchemicalMethod method) {
if (method < 0 || method > 2)
throw OpenMMException("AmoebaVdwForce: Illegal value for alchemical method");
alchemicalMethod = method;
}
void AmoebaVdwForce::setSoftcorePower(double power) {
n = power;
}
double AmoebaVdwForce::getSoftcorePower() const {
return n;
}
void AmoebaVdwForce::setSoftcoreAlpha(double a) {
alpha = a;
}
double AmoebaVdwForce::getSoftcoreAlpha() const {
return alpha;
}
ForceImpl* AmoebaVdwForce::createImpl() const {
return new AmoebaVdwForceImpl(*this);
}
......
plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp
View file @ 36ad9ee7
......@@ -90,10 +90,11 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
map<pair<double, double>, int> classCounts;
for (int i = 0; i < force.getNumParticles(); i++) {
double sigma, epsilon, reduction;
bool isAlchemical;
// The variables reduction, ivindex are not used.
int ivindex;
// Get the sigma and epsilon parameters, ignoring everything else.
force.getParticleParameters(i, ivindex, sigma, epsilon, reduction);
force.getParticleParameters(i, ivindex, sigma, epsilon, reduction, isAlchemical);
pair<double, double> key = make_pair(sigma, epsilon);
map<pair<double, double>, int>::iterator entry = classCounts.find(key);
if (entry == classCounts.end())
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 36ad9ee7
......@@ -2343,9 +2343,10 @@ public:
bool areParticlesIdentical(int particle1, int particle2) {
int iv1, iv2;
double sigma1, sigma2, epsilon1, epsilon2, reduction1, reduction2;
force.getParticleParameters(particle1, iv1, sigma1, epsilon1, reduction1);
force.getParticleParameters(particle2, iv2, sigma2, epsilon2, reduction2);
return (sigma1 == sigma2 && epsilon1 == epsilon2 && reduction1 == reduction2);
bool isAlchemical1, isAlchemical2;
force.getParticleParameters(particle1, iv1, sigma1, epsilon1, reduction1, isAlchemical1);
force.getParticleParameters(particle2, iv2, sigma2, epsilon2, reduction2, isAlchemical2);
return (sigma1 == sigma2 && epsilon1 == epsilon2 && reduction1 == reduction2 && isAlchemical1 == isAlchemical2);
}
private:
const AmoebaVdwForce& force;
......@@ -2378,7 +2379,8 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
for (int i = 0; i < force.getNumParticles(); i++) {
int ivIndex;
double sigma, epsilon, reductionFactor;
force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor);
bool isAlchemical;
force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, isAlchemical);
sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon);
bondReductionAtomsVec[i] = ivIndex;
bondReductionFactorsVec[i] = (float) reductionFactor;
......@@ -2479,7 +2481,8 @@ void CudaCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& context,
for (int i = 0; i < force.getNumParticles(); i++) {
int ivIndex;
double sigma, epsilon, reductionFactor;
force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor);
bool isAlchemical;
force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, isAlchemical);
sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon);
bondReductionAtomsVec[i] = ivIndex;
bondReductionFactorsVec[i] = (float) reductionFactor;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaVdwForce.cpp
View file @ 36ad9ee7
......@@ -140,10 +140,11 @@ void testVdw() {
for (int ii = 0; ii < amoebaVdwForce->getNumParticles(); ii++) {
int indexIV;
double sigma, epsilon, reduction;
amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction);
bool isAlchemical;
amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
sigma *= AngstromToNm;
epsilon *= CalToJoule;
amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction);
amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
}
platformName = "CUDA";
Context context(system, integrator, Platform::getPlatformByName(platformName));
......@@ -170,8 +171,9 @@ void testVdw() {
for (int i = 0; i < numberOfParticles; i++) {
int indexIV;
double mass, sigma, epsilon, reduction;
amoebaVdwForce->getParticleParameters(i, indexIV, sigma, epsilon, reduction);
amoebaVdwForce->setParticleParameters(i, indexIV, 0.9*sigma, 2.0*epsilon, 0.95*reduction);
bool isAlchemical;
amoebaVdwForce->getParticleParameters(i, indexIV, sigma, epsilon, reduction, isAlchemical);
amoebaVdwForce->setParticleParameters(i, indexIV, 0.9*sigma, 2.0*epsilon, 0.95*reduction, isAlchemical);
}
LangevinIntegrator integrator2(0.0, 0.1, 0.01);
Context context2(system, integrator2, Platform::getPlatformByName(platformName));
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ 36ad9ee7
......@@ -996,6 +996,9 @@ ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(std::string
CalcAmoebaVdwForceKernel(name, platform), system(system) {
useCutoff = 0;
usePBC = 0;
alchemicalMethod = 0;
n = 5;
alpha = 0.7;
cutoff = 1.0e+10;
neighborList = NULL;
}
......@@ -1017,14 +1020,16 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
sigmas.resize(numParticles);
epsilons.resize(numParticles);
reductions.resize(numParticles);
isAlchemical.resize(numParticles);
for (int ii = 0; ii < numParticles; ii++) {
int indexIV;
double sigma, epsilon, reduction;
bool alchemical;
std::vector<int> exclusions;
force.getParticleParameters(ii, indexIV, sigma, epsilon, reduction);
force.getParticleParameters(ii, indexIV, sigma, epsilon, reduction, alchemical);
force.getParticleExclusions(ii, exclusions);
for (unsigned int jj = 0; jj < exclusions.size(); jj++) {
allExclusions[ii].insert(exclusions[jj]);
......@@ -1034,6 +1039,7 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
sigmas[ii] = sigma;
epsilons[ii] = epsilon;
reductions[ii] = reduction;
isAlchemical[ii] = alchemical;
}
sigmaCombiningRule = force.getSigmaCombiningRule();
epsilonCombiningRule = force.getEpsilonCombiningRule();
......@@ -1042,7 +1048,9 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
cutoff = force.getCutoffDistance();
neighborList = useCutoff ? new NeighborList() : NULL;
dispersionCoefficient = force.getUseDispersionCorrection() ? AmoebaVdwForceImpl::calcDispersionCorrection(system, force) : 0.0;
alchemicalMethod = (int) force.getAlchemicalMethod();
n = force.getSoftcorePower();
alpha = force.getSoftcoreAlpha();
}
double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
......@@ -1051,6 +1059,7 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
vector<Vec3>& forceData = extractForces(context);
AmoebaReferenceVdwForce vdwForce(sigmaCombiningRule, epsilonCombiningRule);
double energy;
double lambda = context.getParameter("AmoebaVdwLambda");
if (useCutoff) {
vdwForce.setCutoff(cutoff);
computeNeighborListVoxelHash(*neighborList, numParticles, posData, allExclusions, extractBoxVectors(context), usePBC, cutoff, 0.0);
......@@ -1062,14 +1071,16 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
throw OpenMMException("The periodic box size has decreased to less than twice the cutoff.");
}
vdwForce.setPeriodicBox(boxVectors);
energy = vdwForce.calculateForceAndEnergy(numParticles, posData, indexIVs, sigmas, epsilons, reductions, *neighborList, forceData);
energy = vdwForce.calculateForceAndEnergy(numParticles, lambda, posData, indexIVs, sigmas, epsilons,
reductions, isAlchemical, *neighborList, forceData);
energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
} else {
vdwForce.setNonbondedMethod(AmoebaReferenceVdwForce::CutoffNonPeriodic);
}
} else {
vdwForce.setNonbondedMethod(AmoebaReferenceVdwForce::NoCutoff);
energy = vdwForce.calculateForceAndEnergy(numParticles, posData, indexIVs, sigmas, epsilons, reductions, allExclusions, forceData);
energy = vdwForce.calculateForceAndEnergy(numParticles, lambda, posData, indexIVs, sigmas, epsilons,
reductions, isAlchemical, allExclusions, forceData);
}
return static_cast<double>(energy);
}
......@@ -1083,11 +1094,13 @@ void ReferenceCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& con
for (int i = 0; i < numParticles; ++i) {
int indexIV;
double sigma, epsilon, reduction;
force.getParticleParameters(i, indexIV, sigma, epsilon, reduction);
bool alchemical;
force.getParticleParameters(i, indexIV, sigma, epsilon, reduction, alchemical);
indexIVs[i] = indexIV;
sigmas[i] = sigma;
epsilons[i] = epsilon;
reductions[i]= reduction;
isAlchemical[i]= alchemical;
}
}
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ 36ad9ee7
......@@ -499,6 +499,9 @@ private:
int numParticles;
int useCutoff;
int usePBC;
int alchemicalMethod;
double n;
double alpha;
double cutoff;
double dispersionCoefficient;
std::vector<int> indexIVs;
......@@ -506,6 +509,7 @@ private:
std::vector<double> sigmas;
std::vector<double> epsilons;
std::vector<double> reductions;
std::vector<bool> isAlchemical;
std::string sigmaCombiningRule;
std::string epsilonCombiningRule;
const System& system;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp
View file @ 36ad9ee7
......@@ -32,7 +32,7 @@
using std::vector;
using namespace OpenMM;
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce() : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9) {
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce() : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9), _n(5), _alpha(0.7), _alchemicalMethod(None) {
setTaperCoefficients(_cutoff);
setSigmaCombiningRule("ARITHMETIC");
......@@ -40,7 +40,7 @@ AmoebaReferenceVdwForce::AmoebaReferenceVdwForce() : _nonbondedMethod(NoCutoff),
}
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule) : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9) {
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule) : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9), _n(5), _alpha(0.7), _alchemicalMethod(None) {
setTaperCoefficients(_cutoff);
setSigmaCombiningRule(sigmaCombiningRule);
......@@ -55,6 +55,32 @@ void AmoebaReferenceVdwForce::setNonbondedMethod(AmoebaReferenceVdwForce::Nonbon
_nonbondedMethod = nonbondedMethod;
}
AmoebaReferenceVdwForce::AlchemicalMethod AmoebaReferenceVdwForce::getAlchemicalMethod() const {
return _alchemicalMethod;
}
void AmoebaReferenceVdwForce::setAlchemicalMethod(AmoebaReferenceVdwForce::AlchemicalMethod alchemicalMethod){
_alchemicalMethod = alchemicalMethod;
}
void AmoebaReferenceVdwForce::setSoftcorePower(int n) {
_n = n;
}
int AmoebaReferenceVdwForce::getSoftcorePower() const {
return _n;
}
void AmoebaReferenceVdwForce::setSoftcoreAlpha(double alpha) {
_alpha = alpha;
}
double AmoebaReferenceVdwForce::getSoftcoreAlpha() const {
return _alpha;
}
void AmoebaReferenceVdwForce::setTaperCoefficients(double cutoff) {
_taperCutoff = cutoff*_taperCutoffFactor;
if (_taperCutoff != cutoff) {
......@@ -170,13 +196,15 @@ void AmoebaReferenceVdwForce::addReducedForce(unsigned int particleI, unsigned i
forces[particleIV][2] += sign*force[2]*(1.0 - reduction);
}
double AmoebaReferenceVdwForce::calculatePairIxn(double combinedSigma, double combinedEpsilon,
double AmoebaReferenceVdwForce::calculatePairIxn(double combinedSigma, double combinedEpsilon, double softcore,
const Vec3& particleIPosition,
const Vec3& particleJPosition,
Vec3& force) const {
static const double dhal = 0.07;
static const double ghal = 0.12;
static const double dhal1 = 1.07;
static const double ghal1 = 1.12;
// get deltaR, R2, and R between 2 atoms
......@@ -188,29 +216,25 @@ double AmoebaReferenceVdwForce::calculatePairIxn(double combinedSigma, double co
double r_ij_2 = deltaR[ReferenceForce::R2Index];
double r_ij = deltaR[ReferenceForce::RIndex];
double sigma_7 = combinedSigma*combinedSigma*combinedSigma;
sigma_7 = sigma_7*sigma_7*combinedSigma;
double r_ij_6 = r_ij_2*r_ij_2*r_ij_2;
double r_ij_7 = r_ij_6*r_ij;
double rho = r_ij_7 + ghal*sigma_7;
double tau = (dhal + 1.0)/(r_ij + dhal*combinedSigma);
double tau_7 = tau*tau*tau;
tau_7 = tau_7*tau_7*tau;
double dtau = tau/(dhal + 1.0);
double ratio = (sigma_7/rho);
double gtau = combinedEpsilon*tau_7*r_ij_6*(ghal+1.0)*ratio*ratio;
double energy = combinedEpsilon*tau_7*sigma_7*((ghal+1.0)*sigma_7/rho - 2.0);
double dEdR = -7.0*(dtau*energy + gtau);
double rho = r_ij / combinedSigma;
double rho2 = rho * rho;
double rho6 = rho2 * rho2 * rho2;
double rhoplus = rho + dhal;
double rhodec2 = rhoplus * rhoplus;
double rhodec = rhodec2 * rhodec2 * rhodec2;
double s1 = 1.0 / (softcore + rhodec * rhoplus);
double s2 = 1.0 / (softcore + rho6 * rho + 0.12);
double point72 = dhal1 * dhal1;
double t1 = dhal1 * point72 * point72 * point72 * s1;
double t2 = ghal1 * s2;
double t2min = t2 - 2;
double dt1 = -7.0 * rhodec * t1 * s1;
double dt2 = -7.0 * rho6 * t2 * s2;
double energy = combinedEpsilon * t1 * t2min;
double dEdR = combinedEpsilon * (dt1 * t2min + t1 * dt2) / combinedSigma;
// tapering
if ((_nonbondedMethod == CutoffNonPeriodic || _nonbondedMethod == CutoffPeriodic) && r_ij > _taperCutoff) {
double delta = r_ij - _taperCutoff;
double taper = 1.0 + delta*delta*delta*(_taperCoefficients[C3] + delta*(_taperCoefficients[C4] + delta*_taperCoefficients[C5]));
......@@ -248,12 +272,13 @@ void AmoebaReferenceVdwForce::setReducedPositions(int numParticles,
}
}
double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double lambda,
const vector<Vec3>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<double>& sigmas,
const std::vector<double>& epsilons,
const std::vector<double>& reductions,
const std::vector<bool>& isAlchemical,
const std::vector< std::set<int> >& allExclusions,
vector<Vec3>& forces) const {
......@@ -277,6 +302,7 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
double sigmaI = sigmas[ii];
double epsilonI = epsilons[ii];
bool isAlchemicalI = isAlchemical[ii];
for (int jj : allExclusions[ii])
exclusions[jj] = 1;
......@@ -285,9 +311,23 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
double combinedSigma = (this->*_combineSigmas)(sigmaI, sigmas[jj]);
double combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj]);
double softcore = 0.0;
// Apply softcore modifications using an alchemical method.
if (this->_alchemicalMethod != None) {
if (isAlchemicalI || isAlchemical[jj]) {
// Decouple maintains interactions between two softcore atoms.
if (this->_alchemicalMethod == Decouple && isAlchemicalI && isAlchemical[jj]) {
// No softcore
} else {
combinedEpsilon *= pow(lambda, this->_n);
softcore = this->_alpha * pow(1.0 - lambda, 2);
}
}
}
Vec3 force;
energy += calculatePairIxn(combinedSigma, combinedEpsilon,
energy += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
reducedPositions[ii], reducedPositions[jj],
force);
......@@ -316,12 +356,13 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
return energy;
}
double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double lambda,
const vector<Vec3>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<double>& sigmas,
const std::vector<double>& epsilons,
const std::vector<double>& reductions,
const std::vector<bool>& isAlchemical,
const NeighborList& neighborList,
vector<Vec3>& forces) const {
......@@ -345,9 +386,25 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
double combinedSigma = (this->*_combineSigmas)(sigmas[siteI], sigmas[siteJ]);
double combinedEpsilon = (this->*_combineEpsilons)(epsilons[siteI], epsilons[siteJ]);
double softcore = 0.0;
int isAlchemicalI = isAlchemical[siteI];
int isAlchemicalJ = isAlchemical[siteJ];
// Apply softcore modifications using an alchemical method.
if (this->_alchemicalMethod != None) {
if (isAlchemicalI || isAlchemicalJ) {
// Decouple maintains interactions between two softcore atoms.
if (this->_alchemicalMethod == Decouple && isAlchemicalI && isAlchemicalJ) {
// No softcore
} else {
combinedEpsilon *= pow(lambda, this->_n);
softcore = this->_alpha * pow(1.0 - lambda, 2);
}
}
}
Vec3 force;
energy += calculatePairIxn(combinedSigma, combinedEpsilon,
energy += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
reducedPositions[siteI], reducedPositions[siteJ], force);
if (indexIVs[siteI] == siteI) {
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.h
View file @ 36ad9ee7
......@@ -63,6 +63,25 @@ public:
*/
CutoffPeriodic = 2,
};
/**
* This is an enumeration of the different alchemical methods used when applying softcore interactions.
*/
enum AlchemicalMethod {
/**
* All vdW interactions are treated normally. This is the default.
*/
None = 0,
/**
* Maintain full strength vdW interactions between two alchemical particles.
*/
Decouple = 1,
/**
* Interactions between two alchemical particles are turned off at lambda=0.
*/
Annihilate = 2,
};
/**---------------------------------------------------------------------------------------
......@@ -109,6 +128,27 @@ public:
void setNonbondedMethod(NonbondedMethod nonbondedMethod);
/**---------------------------------------------------------------------------------------
Get alchemical method
@return alchemical method
--------------------------------------------------------------------------------------- */
AlchemicalMethod getAlchemicalMethod() const;
/**---------------------------------------------------------------------------------------
Set alchemical method
@param alchemical method
--------------------------------------------------------------------------------------- */
void setAlchemicalMethod(AlchemicalMethod alchemicalMethod);
/**---------------------------------------------------------------------------------------
Get cutoff
......@@ -129,6 +169,47 @@ public:
void setCutoff(double cutoff);
/**---------------------------------------------------------------------------------------
Get softcore power
@return n
--------------------------------------------------------------------------------------- */
int getSoftcorePower() const;
/**---------------------------------------------------------------------------------------
Set softcore power
@param n
--------------------------------------------------------------------------------------- */
void setSoftcorePower(int n);
/**---------------------------------------------------------------------------------------
Get softcore alpha
@return alpha
--------------------------------------------------------------------------------------- */
double getSoftcoreAlpha() const;
/**---------------------------------------------------------------------------------------
Set softcore alpha
@param alpha
--------------------------------------------------------------------------------------- */
void setSoftcoreAlpha(double alpha);
/**---------------------------------------------------------------------------------------
Set sigma combining rule
......@@ -184,11 +265,13 @@ public:
Calculate Amoeba Hal vdw ixns
@param numParticles number of particles
@param lambda lambda value
@param particlePositions Cartesian coordinates of particles
@param indexIVs position index for associated reducing particle
@param sigmas particle sigmas
@param epsilons particle epsilons
@param reductions particle reduction factors
@param isAlchemical particle alchemical flag
@param vdwExclusions particle exclusions
@param forces add forces to this vector
......@@ -196,10 +279,11 @@ public:
--------------------------------------------------------------------------------------- */
double calculateForceAndEnergy(int numParticles, const std::vector<OpenMM::Vec3>& particlePositions,
double calculateForceAndEnergy(int numParticles, double lambda, const std::vector<OpenMM::Vec3>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<double>& sigmas, const std::vector<double>& epsilons,
const std::vector<double>& reductions,
const std::vector<bool>& isAlchemical,
const std::vector< std::set<int> >& vdwExclusions,
std::vector<OpenMM::Vec3>& forces) const;
......@@ -208,11 +292,13 @@ public:
Calculate Vdw ixn using neighbor list
@param numParticles number of particles
@param lambda lambda value
@param particlePositions Cartesian coordinates of particles
@param indexIVs position index for associated reducing particle
@param sigmas particle sigmas
@param epsilons particle epsilons
@param reductions particle reduction factors
@param isAlchemical particle alchemical flag
@param neighborList neighbor list
@param forces add forces to this vector
......@@ -220,17 +306,16 @@ public:
--------------------------------------------------------------------------------------- */
double calculateForceAndEnergy(int numParticles, const std::vector<OpenMM::Vec3>& particlePositions,
double calculateForceAndEnergy(int numParticles, double lambda, const std::vector<OpenMM::Vec3>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<double>& sigmas, const std::vector<double>& epsilons,
const std::vector<double>& reductions,
const std::vector<bool>& isAlchemical,
const NeighborList& neighborList,
std::vector<OpenMM::Vec3>& forces) const;
private:
// taper coefficient indices
static const int C3=0;
static const int C4=1;
static const int C5=2;
......@@ -238,6 +323,9 @@ private:
std::string _sigmaCombiningRule;
std::string _epsilonCombiningRule;
NonbondedMethod _nonbondedMethod;
AlchemicalMethod _alchemicalMethod;
double _n;
double _alpha;
double _cutoff;
double _taperCutoffFactor;
double _taperCutoff;
......@@ -245,14 +333,13 @@ private:
Vec3 _periodicBoxVectors[3];
CombiningFunction _combineSigmas;
double arithmeticSigmaCombiningRule(double sigmaI, double sigmaJ) const;
double geometricSigmaCombiningRule(double sigmaI, double sigmaJ) const;
double cubicMeanSigmaCombiningRule(double sigmaI, double sigmaJ) const;
double geometricSigmaCombiningRule(double sigmaI, double sigmaJ) const;
double cubicMeanSigmaCombiningRule(double sigmaI, double sigmaJ) const;
CombiningFunction _combineEpsilons;
double arithmeticEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
double geometricEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
double harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
double hhgEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
double geometricEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
double harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
double hhgEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
/**---------------------------------------------------------------------------------------
......@@ -308,6 +395,7 @@ private:
@param combindedSigma combined sigmas
@param combindedEpsilon combined epsilons
@param softcore softcore offset parameter
@param particleIPosition particle I position
@param particleJPosition particle J position
@param force output force
......@@ -316,7 +404,7 @@ private:
--------------------------------------------------------------------------------------- */
double calculatePairIxn(double combindedSigma, double combindedEpsilon,
double calculatePairIxn(double combindedSigma, double combindedEpsilon, double softcore,
const Vec3& particleIPosition, const Vec3& particleJPosition,
Vec3& force) const;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaVdwForce.cpp
View file @ 36ad9ee7
......@@ -142,10 +142,11 @@ void testVdw() {
for (int ii = 0; ii < amoebaVdwForce->getNumParticles(); ii++) {
int indexIV;
double sigma, epsilon, reduction;
amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction);
bool isAlchemical;
amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
sigma *= AngstromToNm;
epsilon *= CalToJoule;
amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction);
amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
}
platformName = "Reference";
Context context(system, integrator, Platform::getPlatformByName(platformName));
......@@ -173,8 +174,9 @@ void testVdw() {
for (int i = 0; i < numberOfParticles; i++) {
int indexIV;
double mass, sigma, epsilon, reduction;
amoebaVdwForce->getParticleParameters(i, indexIV, sigma, epsilon, reduction);
amoebaVdwForce->setParticleParameters(i, indexIV, 0.9*sigma, 2.0*epsilon, 0.95*reduction);
bool isAlchemical;
amoebaVdwForce->getParticleParameters(i, indexIV, sigma, epsilon, reduction, isAlchemical);
amoebaVdwForce->setParticleParameters(i, indexIV, 0.9*sigma, 2.0*epsilon, 0.95*reduction, isAlchemical);
}
LangevinIntegrator integrator2(0.0, 0.1, 0.01);
Context context2(system, integrator2, Platform::getPlatformByName(platformName));
......
plugins/amoeba/serialization/src/AmoebaVdwForceProxy.cpp
View file @ 36ad9ee7
......@@ -49,18 +49,22 @@ void AmoebaVdwForceProxy::serialize(const void* object, SerializationNode& node)
node.setStringProperty("SigmaCombiningRule", force.getSigmaCombiningRule());
node.setStringProperty("EpsilonCombiningRule", force.getEpsilonCombiningRule());
node.setDoubleProperty("VdwCutoff", force.getCutoffDistance());
node.setIntProperty("method", (int) force.getNonbondedMethod());
node.setDoubleProperty("n", force.getSoftcorePower());
node.setDoubleProperty("alpha", force.getSoftcoreAlpha());
node.setIntProperty("alchemicalMethod", (int) force.getAlchemicalMethod());
SerializationNode& particles = node.createChildNode("VdwParticles");
for (unsigned int ii = 0; ii < static_cast<unsigned int>(force.getNumParticles()); ii++) {
int ivIndex;
double sigma, epsilon, reductionFactor;
force.getParticleParameters(ii, ivIndex, sigma, epsilon, reductionFactor);
bool isAlchemical;
force.getParticleParameters(ii, ivIndex, sigma, epsilon, reductionFactor, isAlchemical);
SerializationNode& particle = particles.createChildNode("Particle");
particle.setIntProperty("ivIndex", ivIndex).setDoubleProperty("sigma", sigma).setDoubleProperty("epsilon", epsilon).setDoubleProperty("reductionFactor", reductionFactor);
particle.setIntProperty("ivIndex", ivIndex).setDoubleProperty("sigma", sigma).setDoubleProperty("epsilon", epsilon).setDoubleProperty("reductionFactor", reductionFactor).setBoolProperty("isAlchemical",isAlchemical);
std::vector< int > exclusions;
force.getParticleExclusions(ii, exclusions);
......@@ -85,10 +89,14 @@ void* AmoebaVdwForceProxy::deserialize(const SerializationNode& node) const {
force->setCutoffDistance(node.getDoubleProperty("VdwCutoff"));
force->setNonbondedMethod((AmoebaVdwForce::NonbondedMethod) node.getIntProperty("method"));
force->setAlchemicalMethod((AmoebaVdwForce::AlchemicalMethod) node.getIntProperty("alchemicalMethod"));
force->setSoftcorePower(node.getDoubleProperty("n"));
force->setSoftcoreAlpha(node.getDoubleProperty("alpha"));
const SerializationNode& particles = node.getChildNode("VdwParticles");
for (unsigned int ii = 0; ii < particles.getChildren().size(); ii++) {
const SerializationNode& particle = particles.getChildren()[ii];
force->addParticle(particle.getIntProperty("ivIndex"), particle.getDoubleProperty("sigma"), particle.getDoubleProperty("epsilon"), particle.getDoubleProperty("reductionFactor"));
force->addParticle(particle.getIntProperty("ivIndex"), particle.getDoubleProperty("sigma"), particle.getDoubleProperty("epsilon"), particle.getDoubleProperty("reductionFactor"), particle.getBoolProperty("isAlchemical"));
// exclusions
......
plugins/amoeba/serialization/tests/TestSerializeAmoebaVdwForce.cpp
View file @ 36ad9ee7
......@@ -50,10 +50,11 @@ void testSerialization() {
force1.setEpsilonCombiningRule("GEOMETRIC");
force1.setCutoff(0.9);
force1.setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic);
force1.setAlchemicalMethod(AmoebaVdwForce::None);
force1.addParticle(0, 1.0, 2.0, 0.9);
force1.addParticle(1, 1.1, 2.1, 0.9);
force1.addParticle(2, 1.3, 4.1, 0.9);
force1.addParticle(0, 1.0, 2.0, 0.9, false);
force1.addParticle(1, 1.1, 2.1, 0.9, true);
force1.addParticle(2, 1.3, 4.1, 0.9, false);
for (unsigned int ii = 0; ii < 3; ii++) {
std::vector< int > exclusions;
exclusions.push_back(ii);
......@@ -76,6 +77,7 @@ void testSerialization() {
ASSERT_EQUAL(force1.getEpsilonCombiningRule(), force2.getEpsilonCombiningRule());
ASSERT_EQUAL(force1.getCutoff(), force2.getCutoff());
ASSERT_EQUAL(force1.getNonbondedMethod(), force2.getNonbondedMethod());
ASSERT_EQUAL(force1.getAlchemicalMethod(), force2.getAlchemicalMethod());
ASSERT_EQUAL(force1.getNumParticles(), force2.getNumParticles());
......@@ -87,13 +89,17 @@ void testSerialization() {
double sigma1, epsilon1, reductionFactor1;
double sigma2, epsilon2, reductionFactor2;
force1.getParticleParameters(ii, ivIndex1, sigma1, epsilon1, reductionFactor1);
force2.getParticleParameters(ii, ivIndex2, sigma2, epsilon2, reductionFactor2);
bool isAlchemical1;
bool isAlchemical2;
force1.getParticleParameters(ii, ivIndex1, sigma1, epsilon1, reductionFactor1, isAlchemical1);
force2.getParticleParameters(ii, ivIndex2, sigma2, epsilon2, reductionFactor2, isAlchemical2);
ASSERT_EQUAL(ivIndex1, ivIndex2);
ASSERT_EQUAL(sigma1, sigma2);
ASSERT_EQUAL(epsilon1, epsilon2);
ASSERT_EQUAL(reductionFactor1, reductionFactor2);
ASSERT_EQUAL(isAlchemical1, isAlchemical2);
}
for (unsigned int ii = 0; ii < static_cast<unsigned int>(force1.getNumParticles()); ii++) {
......
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