Skip to content

GitLab

Unverified Commit 68424035 authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #2377 from mjschnie/master 

Add Softcore support to AmoebaVdwForce.
parents 5bc1013f 926e7b9a
Hide whitespace changes
Inline Side-by-side
Showing with 665 additions and 132 deletions
plugins/amoeba/openmmapi/include/openmm/AmoebaVdwForce.h
View file @ 68424035
...@@ -50,10 +50,31 @@ namespace OpenMM { ...@@ -50,10 +50,31 @@ namespace OpenMM {
* particle to another one. This is typically done for hydrogens to place the interaction site slightly * particle to another one. This is typically done for hydrogens to place the interaction site slightly
* closer to the parent atom. The fraction is known as the "reduction factor", since it reduces the distance * closer to the parent atom. The fraction is known as the "reduction factor", since it reduces the distance
* from the parent atom to the interaction site. * from the parent atom to the interaction site.
*
* Support is also available for softcore interactions based on setting a per particle alchemical flag and
* setting the AmoebaVdwForce to use an "AlchemicalMethod" -- either Decouple or Annihilate.
* For Decouple, two alchemical atoms interact normally. For Annihilate, all interactions involving an
* alchemical atom are influenced. The softcore state is specified by setting a single
* Context parameter "AmoebaVdwLambda" between 0.0 and 1.0.
*
* The softcore functional form can be modified by setting the softcore power (default of 5) and the softcore
* alpha (default of 0,7). For more information on the softcore functional form see Eq. 2 from:
* Jiao, D.; Golubkov, P. A.; Darden, T. A.; Ren, P.,
* Calculation of protein-ligand binding free energy by using a polarizable potential.
* Proc. Natl. Acad. Sci. U.S.A. 2008, 105 (17), 6290-6295.
* https://www.pnas.org/content/105/17/6290.
*/ */
class OPENMM_EXPORT_AMOEBA AmoebaVdwForce : public Force { class OPENMM_EXPORT_AMOEBA AmoebaVdwForce : public Force {
public: public:
/**
* This is the name of the parameter which stores the current Amoeba vdW lambda value.
*/
static const std::string& Lambda() {
static const std::string key = "AmoebaVdwLambda";
return key;
}
/** /**
* This is an enumeration of the different methods that may be used for handling long range nonbonded forces. * This is an enumeration of the different methods that may be used for handling long range nonbonded forces.
*/ */
...@@ -70,6 +91,24 @@ public: ...@@ -70,6 +91,24 @@ public:
CutoffPeriodic = 1, CutoffPeriodic = 1,
}; };
/**
* 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,
};
/** /**
* Create an Amoeba VdwForce. * Create an Amoeba VdwForce.
*/ */
...@@ -91,8 +130,10 @@ public: ...@@ -91,8 +130,10 @@ public:
* @param epsilon vdw epsilon * @param epsilon vdw epsilon
* @param reductionFactor the fraction of the distance along the line from the parent particle to this particle * @param reductionFactor the fraction of the distance along the line from the parent particle to this particle
* at which the interaction site should be placed * at which the interaction site should be placed
* @param isAlchemical if true, this vdW particle is undergoing an alchemical change.
*/ */
void setParticleParameters(int particleIndex, int parentIndex, double sigma, double epsilon, double reductionFactor); void setParticleParameters(int particleIndex, int parentIndex, double sigma, double epsilon,
double reductionFactor, bool isAlchemical = false);
/** /**
* Get the force field parameters for a vdw particle. * Get the force field parameters for a vdw particle.
...@@ -103,8 +144,10 @@ public: ...@@ -103,8 +144,10 @@ public:
* @param[out] epsilon vdw epsilon * @param[out] epsilon vdw epsilon
* @param[out] reductionFactor the fraction of the distance along the line from the parent particle to this particle * @param[out] reductionFactor the fraction of the distance along the line from the parent particle to this particle
* at which the interaction site should be placed * at which the interaction site should be placed
* @param[out] isAlchemical if true, this vdW particle is undergoing an alchemical change.
*/ */
void getParticleParameters(int particleIndex, int& parentIndex, double& sigma, double& epsilon, double& reductionFactor) const; void getParticleParameters(int particleIndex, int& parentIndex, double& sigma, double& epsilon,
double& reductionFactor, bool& isAlchemical) const;
/** /**
...@@ -115,9 +158,10 @@ public: ...@@ -115,9 +158,10 @@ public:
* @param epsilon vdw epsilon * @param epsilon vdw epsilon
* @param reductionFactor the fraction of the distance along the line from the parent particle to this particle * @param reductionFactor the fraction of the distance along the line from the parent particle to this particle
* at which the interaction site should be placed * at which the interaction site should be placed
* @param isAlchemical if true, this vdW particle is undergoing an alchemical change.
* @return index of added particle * @return index of added particle
*/ */
int addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor); int addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor, bool isAlchemical = false);
/** /**
* Set sigma combining rule * Set sigma combining rule
...@@ -223,6 +267,37 @@ public: ...@@ -223,6 +267,37 @@ public:
* Set the method used for handling long range nonbonded interactions. * Set the method used for handling long range nonbonded interactions.
*/ */
void setNonbondedMethod(NonbondedMethod method); void setNonbondedMethod(NonbondedMethod method);
/**
* Set the softcore power on lambda (default = 5).
*/
void setSoftcorePower(int n);
/**
* Get the softcore power on lambda.
*/
int getSoftcorePower() const;
/**
* Set the softcore alpha value (default = 0.7).
*/
void setSoftcoreAlpha(double alpha);
/**
* Get the softcore alpha value.
*/
double getSoftcoreAlpha() const;
/**
* Get the method used for alchemical interactions.
*/
AlchemicalMethod getAlchemicalMethod() const;
/**
* Set the method used for handling long range nonbonded interactions.
*/
void setAlchemicalMethod(AlchemicalMethod method);
/** /**
* Update the per-particle parameters in a Context to match those stored in this Force object. This method provides * Update the per-particle parameters in a Context to match those stored in this Force object. This method provides
* an efficient method to update certain parameters in an existing Context without needing to reinitialize it. * an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
...@@ -250,6 +325,9 @@ private: ...@@ -250,6 +325,9 @@ private:
NonbondedMethod nonbondedMethod; NonbondedMethod nonbondedMethod;
double cutoff; double cutoff;
bool useDispersionCorrection; bool useDispersionCorrection;
AlchemicalMethod alchemicalMethod;
int n;
double alpha;
std::string sigmaCombiningRule; std::string sigmaCombiningRule;
std::string epsilonCombiningRule; std::string epsilonCombiningRule;
...@@ -267,14 +345,16 @@ class AmoebaVdwForce::VdwInfo { ...@@ -267,14 +345,16 @@ class AmoebaVdwForce::VdwInfo {
public: public:
int parentIndex; int parentIndex;
double reductionFactor, sigma, epsilon, cutoff; double reductionFactor, sigma, epsilon, cutoff;
bool isAlchemical;
VdwInfo() { VdwInfo() {
parentIndex = -1; parentIndex = -1;
reductionFactor = 0.0; reductionFactor = 0.0;
sigma = 1.0; sigma = 1.0;
epsilon = 0.0; epsilon = 0.0;
isAlchemical = false;
} }
VdwInfo(int parentIndex, double sigma, double epsilon, double reductionFactor) : VdwInfo(int parentIndex, double sigma, double epsilon, double reductionFactor, bool isAlchemical) :
parentIndex(parentIndex), reductionFactor(reductionFactor), sigma(sigma), epsilon(epsilon) { parentIndex(parentIndex), reductionFactor(reductionFactor), sigma(sigma), epsilon(epsilon), isAlchemical(isAlchemical) {
} }
}; };
......
plugins/amoeba/openmmapi/include/openmm/internal/AmoebaVdwForceImpl.h
View file @ 68424035
...@@ -60,7 +60,9 @@ public: ...@@ -60,7 +60,9 @@ public:
} }
double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups); double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
std::map<std::string, double> getDefaultParameters() { std::map<std::string, double> getDefaultParameters() {
return std::map<std::string, double>(); // This force field doesn't define any parameters. std::map<std::string, double> parameters;
parameters[AmoebaVdwForce::Lambda()] = 1.0;
return parameters;
} }
std::vector<std::string> getKernelNames(); std::vector<std::string> getKernelNames();
/** /**
......
plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp
View file @ 68424035
...@@ -38,28 +38,30 @@ using namespace OpenMM; ...@@ -38,28 +38,30 @@ using namespace OpenMM;
using std::string; using std::string;
using std::vector; 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) { int AmoebaVdwForce::addParticle(int parentIndex, double sigma, double epsilon, double reductionFactor, bool isAlchemical) {
parameters.push_back(VdwInfo(parentIndex, sigma, epsilon, reductionFactor)); parameters.push_back(VdwInfo(parentIndex, sigma, epsilon, reductionFactor, isAlchemical));
return parameters.size()-1; return parameters.size()-1;
} }
void AmoebaVdwForce::getParticleParameters(int particleIndex, int& parentIndex, 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; parentIndex = parameters[particleIndex].parentIndex;
sigma = parameters[particleIndex].sigma; sigma = parameters[particleIndex].sigma;
epsilon = parameters[particleIndex].epsilon; epsilon = parameters[particleIndex].epsilon;
reductionFactor = parameters[particleIndex].reductionFactor; reductionFactor = parameters[particleIndex].reductionFactor;
isAlchemical = parameters[particleIndex].isAlchemical;
} }
void AmoebaVdwForce::setParticleParameters(int particleIndex, int parentIndex, 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].parentIndex = parentIndex;
parameters[particleIndex].sigma = sigma; parameters[particleIndex].sigma = sigma;
parameters[particleIndex].epsilon = epsilon; parameters[particleIndex].epsilon = epsilon;
parameters[particleIndex].reductionFactor = reductionFactor; parameters[particleIndex].reductionFactor = reductionFactor;
parameters[particleIndex].isAlchemical = isAlchemical;
} }
void AmoebaVdwForce::setSigmaCombiningRule(const std::string& inputSigmaCombiningRule) { void AmoebaVdwForce::setSigmaCombiningRule(const std::string& inputSigmaCombiningRule) {
...@@ -128,6 +130,33 @@ void AmoebaVdwForce::setNonbondedMethod(NonbondedMethod method) { ...@@ -128,6 +130,33 @@ void AmoebaVdwForce::setNonbondedMethod(NonbondedMethod method) {
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(int power) {
n = power;
}
int AmoebaVdwForce::getSoftcorePower() const {
return n;
}
void AmoebaVdwForce::setSoftcoreAlpha(double a) {
alpha = a;
}
double AmoebaVdwForce::getSoftcoreAlpha() const {
return alpha;
}
ForceImpl* AmoebaVdwForce::createImpl() const { ForceImpl* AmoebaVdwForce::createImpl() const {
return new AmoebaVdwForceImpl(*this); return new AmoebaVdwForceImpl(*this);
} }
......
plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp
View file @ 68424035
...@@ -90,10 +90,11 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const ...@@ -90,10 +90,11 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
map<pair<double, double>, int> classCounts; map<pair<double, double>, int> classCounts;
for (int i = 0; i < force.getNumParticles(); i++) { for (int i = 0; i < force.getNumParticles(); i++) {
double sigma, epsilon, reduction; double sigma, epsilon, reduction;
// The variables reduction, ivindex are not used. bool isAlchemical;
// The variables reduction, ivindex and isAlchemical are not used.
int ivindex; int ivindex;
// Get the sigma and epsilon parameters, ignoring everything else. // 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); pair<double, double> key = make_pair(sigma, epsilon);
map<pair<double, double>, int>::iterator entry = classCounts.find(key); map<pair<double, double>, int>::iterator entry = classCounts.find(key);
if (entry == classCounts.end()) if (entry == classCounts.end())
...@@ -151,8 +152,11 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const ...@@ -151,8 +152,11 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
int ndelta = int(double(nstep) * (range - cut)); int ndelta = int(double(nstep) * (range - cut));
double rdelta = (range - cut) / double(ndelta); double rdelta = (range - cut) / double(ndelta);
double offset = cut - 0.5 * rdelta; double offset = cut - 0.5 * rdelta;
double dhal = 0.07; // This magic number also appears in kCalculateAmoebaCudaVdw14_7.cu
double ghal = 0.12; // This magic number also appears in kCalculateAmoebaCudaVdw14_7.cu // Buffered-14-7 buffering constants
double dhal = 0.07;
double ghal = 0.12;
double elrc = 0.0; // This number is incremented and passed out at the end double elrc = 0.0; // This number is incremented and passed out at the end
double e = 0.0; double e = 0.0;
double sigma, epsilon; // The pairwise sigma and epsilon parameters. double sigma, epsilon; // The pairwise sigma and epsilon parameters.
...@@ -259,3 +263,5 @@ std::vector<std::string> AmoebaVdwForceImpl::getKernelNames() { ...@@ -259,3 +263,5 @@ std::vector<std::string> AmoebaVdwForceImpl::getKernelNames() {
void AmoebaVdwForceImpl::updateParametersInContext(ContextImpl& context) { void AmoebaVdwForceImpl::updateParametersInContext(ContextImpl& context) {
kernel.getAs<CalcAmoebaVdwForceKernel>().copyParametersToContext(context, owner); kernel.getAs<CalcAmoebaVdwForceKernel>().copyParametersToContext(context, owner);
} }
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 68424035
...@@ -2343,43 +2343,57 @@ public: ...@@ -2343,43 +2343,57 @@ public:
bool areParticlesIdentical(int particle1, int particle2) { bool areParticlesIdentical(int particle1, int particle2) {
int iv1, iv2; int iv1, iv2;
double sigma1, sigma2, epsilon1, epsilon2, reduction1, reduction2; double sigma1, sigma2, epsilon1, epsilon2, reduction1, reduction2;
force.getParticleParameters(particle1, iv1, sigma1, epsilon1, reduction1); bool isAlchemical1, isAlchemical2;
force.getParticleParameters(particle2, iv2, sigma2, epsilon2, reduction2); force.getParticleParameters(particle1, iv1, sigma1, epsilon1, reduction1, isAlchemical1);
return (sigma1 == sigma2 && epsilon1 == epsilon2 && reduction1 == reduction2); force.getParticleParameters(particle2, iv2, sigma2, epsilon2, reduction2, isAlchemical2);
return (sigma1 == sigma2 && epsilon1 == epsilon2 && reduction1 == reduction2 && isAlchemical1 == isAlchemical2);
} }
private: private:
const AmoebaVdwForce& force; const AmoebaVdwForce& force;
}; };
CudaCalcAmoebaVdwForceKernel::CudaCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system) : CudaCalcAmoebaVdwForceKernel::CudaCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system) :
CalcAmoebaVdwForceKernel(name, platform), cu(cu), system(system), hasInitializedNonbonded(false), nonbonded(NULL) { CalcAmoebaVdwForceKernel(name, platform), cu(cu), system(system), hasInitializedNonbonded(false), nonbonded(NULL), vdwLambdaPinnedBuffer(NULL) {
} }
CudaCalcAmoebaVdwForceKernel::~CudaCalcAmoebaVdwForceKernel() { CudaCalcAmoebaVdwForceKernel::~CudaCalcAmoebaVdwForceKernel() {
cu.setAsCurrent(); cu.setAsCurrent();
if (nonbonded != NULL) if (nonbonded != NULL)
delete nonbonded; delete nonbonded;
if (vdwLambdaPinnedBuffer != NULL)
cuMemFreeHost(vdwLambdaPinnedBuffer);
} }
void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) { void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) {
cu.setAsCurrent(); cu.setAsCurrent();
sigmaEpsilon.initialize<float2>(cu, cu.getPaddedNumAtoms(), "sigmaEpsilon"); sigmaEpsilon.initialize<float2>(cu, cu.getPaddedNumAtoms(), "sigmaEpsilon");
bondReductionAtoms.initialize<int>(cu, cu.getPaddedNumAtoms(), "bondReductionAtoms"); bondReductionAtoms.initialize<int>(cu, cu.getPaddedNumAtoms(), "bondReductionAtoms");
bondReductionFactors.initialize<float>(cu, cu.getPaddedNumAtoms(), "sigmaEpsilon"); bondReductionFactors.initialize<float>(cu, cu.getPaddedNumAtoms(), "bondReductionFactors");
tempPosq.initialize(cu, cu.getPaddedNumAtoms(), cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4), "tempPosq"); tempPosq.initialize(cu, cu.getPaddedNumAtoms(), cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4), "tempPosq");
tempForces.initialize<long long>(cu, 3*cu.getPaddedNumAtoms(), "tempForces"); tempForces.initialize<long long>(cu, 3*cu.getPaddedNumAtoms(), "tempForces");
// Record atom parameters. // Record atom parameters.
vector<float2> sigmaEpsilonVec(cu.getPaddedNumAtoms(), make_float2(0, 1)); vector<float2> sigmaEpsilonVec(cu.getPaddedNumAtoms(), make_float2(0, 1));
vector<float> isAlchemicalVec(cu.getPaddedNumAtoms(), 0);
vector<int> bondReductionAtomsVec(cu.getPaddedNumAtoms(), 0); vector<int> bondReductionAtomsVec(cu.getPaddedNumAtoms(), 0);
vector<float> bondReductionFactorsVec(cu.getPaddedNumAtoms(), 0); vector<float> bondReductionFactorsVec(cu.getPaddedNumAtoms(), 0);
vector<vector<int> > exclusions(cu.getNumAtoms()); vector<vector<int> > exclusions(cu.getNumAtoms());
// Handle Alchemical parameters.
hasAlchemical = force.getAlchemicalMethod() != AmoebaVdwForce::None;
if (hasAlchemical) {
isAlchemical.initialize<float>(cu, cu.getPaddedNumAtoms(), "isAlchemical");
vdwLambda.initialize<float>(cu, 1, "vdwLambda");
CHECK_RESULT(cuMemHostAlloc(&vdwLambdaPinnedBuffer, sizeof(float), 0), "Error allocating vdwLambda pinned buffer");
}
for (int i = 0; i < force.getNumParticles(); i++) { for (int i = 0; i < force.getNumParticles(); i++) {
int ivIndex; int ivIndex;
double sigma, epsilon, reductionFactor; double sigma, epsilon, reductionFactor;
force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor); bool alchemical;
force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical);
sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon); sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon);
isAlchemicalVec[i] = (alchemical) ? 1.0f : 0.0f;
bondReductionAtomsVec[i] = ivIndex; bondReductionAtomsVec[i] = ivIndex;
bondReductionFactorsVec[i] = (float) reductionFactor; bondReductionFactorsVec[i] = (float) reductionFactor;
force.getParticleExclusions(i, exclusions[i]); force.getParticleExclusions(i, exclusions[i]);
...@@ -2392,13 +2406,22 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba ...@@ -2392,13 +2406,22 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
dispersionCoefficient = AmoebaVdwForceImpl::calcDispersionCorrection(system, force); dispersionCoefficient = AmoebaVdwForceImpl::calcDispersionCorrection(system, force);
else else
dispersionCoefficient = 0.0; dispersionCoefficient = 0.0;
// This force is applied based on modified atom positions, where hydrogens have been moved slightly // This force is applied based on modified atom positions, where hydrogens have been moved slightly
// closer to their parent atoms. We therefore create a separate CudaNonbondedUtilities just for // closer to their parent atoms. We therefore create a separate CudaNonbondedUtilities just for
// this force, so it will have its own neighbor list and interaction kernel. // this force, so it will have its own neighbor list and interaction kernel.
nonbonded = new CudaNonbondedUtilities(cu); nonbonded = new CudaNonbondedUtilities(cu);
nonbonded->addParameter(CudaNonbondedUtilities::ParameterInfo("sigmaEpsilon", "float", 2, sizeof(float2), sigmaEpsilon.getDevicePointer())); nonbonded->addParameter(CudaNonbondedUtilities::ParameterInfo("sigmaEpsilon", "float", 2, sizeof(float2), sigmaEpsilon.getDevicePointer()));
if (hasAlchemical) {
isAlchemical.upload(isAlchemicalVec);
((float*) vdwLambdaPinnedBuffer)[0] = 1.0f;
currentVdwLambda = 1.0f;
vdwLambda.upload(vdwLambdaPinnedBuffer, false);
nonbonded->addParameter(CudaNonbondedUtilities::ParameterInfo("isAlchemical", "float", 1, sizeof(float), isAlchemical.getDevicePointer()));
nonbonded->addArgument(CudaNonbondedUtilities::ParameterInfo("vdwLambda", "float", 1, sizeof(float), vdwLambda.getDevicePointer()));
}
// Create the interaction kernel. // Create the interaction kernel.
...@@ -2417,12 +2440,17 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba ...@@ -2417,12 +2440,17 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
replacements["EPSILON_COMBINING_RULE"] = "1"; replacements["EPSILON_COMBINING_RULE"] = "1";
else if (epsilonCombiningRule == "GEOMETRIC") else if (epsilonCombiningRule == "GEOMETRIC")
replacements["EPSILON_COMBINING_RULE"] = "2"; replacements["EPSILON_COMBINING_RULE"] = "2";
else if (epsilonCombiningRule == "HARMONIC") else if (epsilonCombiningRule =="HARMONIC")
replacements["EPSILON_COMBINING_RULE"] = "3"; replacements["EPSILON_COMBINING_RULE"] = "3";
else if (epsilonCombiningRule == "HHG") else if (epsilonCombiningRule == "HHG")
replacements["EPSILON_COMBINING_RULE"] = "4"; replacements["EPSILON_COMBINING_RULE"] = "4";
else else
throw OpenMMException("Illegal combining rule for sigma: "+sigmaCombiningRule); throw OpenMMException("Illegal combining rule for sigma: "+sigmaCombiningRule);
replacements["VDW_ALCHEMICAL_METHOD"] = cu.intToString(force.getAlchemicalMethod());
replacements["VDW_SOFTCORE_POWER"] = cu.intToString(force.getSoftcorePower());
replacements["VDW_SOFTCORE_ALPHA"] = cu.doubleToString(force.getSoftcoreAlpha());
double cutoff = force.getCutoffDistance(); double cutoff = force.getCutoffDistance();
double taperCutoff = cutoff*0.9; double taperCutoff = cutoff*0.9;
replacements["CUTOFF_DISTANCE"] = cu.doubleToString(force.getCutoffDistance()); replacements["CUTOFF_DISTANCE"] = cu.doubleToString(force.getCutoffDistance());
...@@ -2449,6 +2477,17 @@ double CudaCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeF ...@@ -2449,6 +2477,17 @@ double CudaCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeF
hasInitializedNonbonded = true; hasInitializedNonbonded = true;
nonbonded->initialize(system); nonbonded->initialize(system);
} }
if (hasAlchemical) {
float contextLambda = context.getParameter(AmoebaVdwForce::Lambda());
if (contextLambda != currentVdwLambda) {
// Non-blocking copy of vdwLambda to device memory
((float*) vdwLambdaPinnedBuffer)[0] = contextLambda;
vdwLambda.upload(vdwLambdaPinnedBuffer, false);
currentVdwLambda = contextLambda;
}
}
cu.getPosq().copyTo(tempPosq); cu.getPosq().copyTo(tempPosq);
cu.getForce().copyTo(tempForces); cu.getForce().copyTo(tempForces);
void* prepareArgs[] = {&cu.getForce().getDevicePointer(), &cu.getPosq().getDevicePointer(), &tempPosq.getDevicePointer(), void* prepareArgs[] = {&cu.getForce().getDevicePointer(), &cu.getPosq().getDevicePointer(), &tempPosq.getDevicePointer(),
...@@ -2472,19 +2511,22 @@ void CudaCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& context, ...@@ -2472,19 +2511,22 @@ void CudaCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& context,
throw OpenMMException("updateParametersInContext: The number of particles has changed"); throw OpenMMException("updateParametersInContext: The number of particles has changed");
// Record the per-particle parameters. // Record the per-particle parameters.
vector<float2> sigmaEpsilonVec(cu.getPaddedNumAtoms(), make_float2(0, 1)); vector<float2> sigmaEpsilonVec(cu.getPaddedNumAtoms(), make_float2(0, 1));
vector<float> isAlchemicalVec(cu.getPaddedNumAtoms(), 0);
vector<int> bondReductionAtomsVec(cu.getPaddedNumAtoms(), 0); vector<int> bondReductionAtomsVec(cu.getPaddedNumAtoms(), 0);
vector<float> bondReductionFactorsVec(cu.getPaddedNumAtoms(), 0); vector<float> bondReductionFactorsVec(cu.getPaddedNumAtoms(), 0);
for (int i = 0; i < force.getNumParticles(); i++) { for (int i = 0; i < force.getNumParticles(); i++) {
int ivIndex; int ivIndex;
double sigma, epsilon, reductionFactor; double sigma, epsilon, reductionFactor;
force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor); bool alchemical;
force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor, alchemical);
sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon); sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon);
isAlchemicalVec[i] = (alchemical) ? 1.0f : 0.0f;
bondReductionAtomsVec[i] = ivIndex; bondReductionAtomsVec[i] = ivIndex;
bondReductionFactorsVec[i] = (float) reductionFactor; bondReductionFactorsVec[i] = (float) reductionFactor;
} }
sigmaEpsilon.upload(sigmaEpsilonVec); sigmaEpsilon.upload(sigmaEpsilonVec);
if (hasAlchemical) isAlchemical.upload(isAlchemicalVec);
bondReductionAtoms.upload(bondReductionAtomsVec); bondReductionAtoms.upload(bondReductionAtomsVec);
bondReductionFactors.upload(bondReductionFactorsVec); bondReductionFactors.upload(bondReductionFactorsVec);
if (force.getUseDispersionCorrection()) if (force.getUseDispersionCorrection())
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h
View file @ 68424035
...@@ -547,7 +547,7 @@ public: ...@@ -547,7 +547,7 @@ public:
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the AmoebaMultipoleForce this kernel will be used for * @param force the AmoebaVdwForce this kernel will be used for
*/ */
void initialize(const System& system, const AmoebaVdwForce& force); void initialize(const System& system, const AmoebaVdwForce& force);
/** /**
...@@ -571,6 +571,18 @@ private: ...@@ -571,6 +571,18 @@ private:
CudaContext& cu; CudaContext& cu;
const System& system; const System& system;
bool hasInitializedNonbonded; bool hasInitializedNonbonded;
// True if the AmoebaVdwForce AlchemicalMethod is not None.
bool hasAlchemical;
// Pinned host memory; allocated if necessary in initialize, and freed in the destructor.
void* vdwLambdaPinnedBuffer;
// Device memory for the alchemical state.
CudaArray vdwLambda;
// Only update device memory when lambda changes.
float currentVdwLambda;
// Per particle alchemical flag.
CudaArray isAlchemical;
double dispersionCoefficient; double dispersionCoefficient;
CudaArray sigmaEpsilon; CudaArray sigmaEpsilon;
CudaArray bondReductionAtoms; CudaArray bondReductionAtoms;
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaVdwForce2.cu
View file @ 68424035
...@@ -26,20 +26,38 @@ ...@@ -26,20 +26,38 @@
real epsilon_s = SQRT(sigmaEpsilon1.y) + SQRT(sigmaEpsilon2.y); real epsilon_s = SQRT(sigmaEpsilon1.y) + SQRT(sigmaEpsilon2.y);
real epsilon = (epsilon_s == 0.0f ? (real) 0 : 4*sigmaEpsilon1.y*sigmaEpsilon2.y/(epsilon_s*epsilon_s)); real epsilon = (epsilon_s == 0.0f ? (real) 0 : 4*sigmaEpsilon1.y*sigmaEpsilon2.y/(epsilon_s*epsilon_s));
#endif #endif
real r6 = r2*r2*r2; real softcore = 0.0f;
real r7 = r6*r; #if VDW_ALCHEMICAL_METHOD == 1
real sigma7 = sigma*sigma; if (isAlchemical1 != isAlchemical2) {
sigma7 = sigma7*sigma7*sigma7*sigma; #elif VDW_ALCHEMICAL_METHOD == 2
real rho = r7 + sigma7*0.12f; if (isAlchemical1 || isAlchemical2) {
real invRho = RECIP(rho); #endif
real tau = 1.07f/(r + 0.07f*sigma); #if VDW_ALCHEMICAL_METHOD != 0
real tau7 = tau*tau*tau; real lambda = vdwLambda[0];
tau7 = tau7*tau7*tau; epsilon = epsilon * POW(lambda, VDW_SOFTCORE_POWER);
real dTau = tau/1.07f; softcore = VDW_SOFTCORE_ALPHA * (1.0f - lambda) * (1.0f - lambda);
real tmp = sigma7*invRho; }
real gTau = epsilon*tau7*r6*1.12f*tmp*tmp; #endif
real termEnergy = epsilon*sigma7*tau7*((sigma7*1.12f*invRho)-2.0f); real dhal = 0.07f;
real deltaE = -7.0f*(dTau*termEnergy+gTau); real ghal = 0.12f;
real dhal1 = 1.07f;
real ghal1 = 1.12f;
real rho = r / sigma;
real rho2 = rho * rho;
real rho6 = rho2 * rho2 * rho2;
real rhoplus = rho + dhal;
real rhodec2 = rhoplus * rhoplus;
real rhodec = rhodec2 * rhodec2 * rhodec2;
real s1 = 1.0f / (softcore + rhodec * rhoplus);
real s2 = 1.0f / (softcore + rho6 * rho + ghal);
real point72 = dhal1 * dhal1;
real t1 = dhal1 * point72 * point72 * point72 * s1;
real t2 = ghal1 * s2;
real t2min = t2 - 2.0f;
real dt1 = -7.0f * rhodec * t1 * s1;
real dt2 = -7.0f * rho6 * t2 * s2;
real termEnergy = epsilon * t1 * t2min;
real deltaE = epsilon * (dt1 * t2min + t1 * dt2) / sigma;
#ifdef USE_CUTOFF #ifdef USE_CUTOFF
if (r > TAPER_CUTOFF) { if (r > TAPER_CUTOFF) {
real x = r-TAPER_CUTOFF; real x = r-TAPER_CUTOFF;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaVdwForce.cpp
View file @ 68424035
...@@ -140,10 +140,11 @@ void testVdw() { ...@@ -140,10 +140,11 @@ void testVdw() {
for (int ii = 0; ii < amoebaVdwForce->getNumParticles(); ii++) { for (int ii = 0; ii < amoebaVdwForce->getNumParticles(); ii++) {
int indexIV; int indexIV;
double sigma, epsilon, reduction; double sigma, epsilon, reduction;
amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction); bool isAlchemical;
amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
sigma *= AngstromToNm; sigma *= AngstromToNm;
epsilon *= CalToJoule; epsilon *= CalToJoule;
amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction); amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
} }
platformName = "CUDA"; platformName = "CUDA";
Context context(system, integrator, Platform::getPlatformByName(platformName)); Context context(system, integrator, Platform::getPlatformByName(platformName));
...@@ -170,8 +171,9 @@ void testVdw() { ...@@ -170,8 +171,9 @@ void testVdw() {
for (int i = 0; i < numberOfParticles; i++) { for (int i = 0; i < numberOfParticles; i++) {
int indexIV; int indexIV;
double mass, sigma, epsilon, reduction; double mass, sigma, epsilon, reduction;
amoebaVdwForce->getParticleParameters(i, indexIV, sigma, epsilon, reduction); bool isAlchemical;
amoebaVdwForce->setParticleParameters(i, indexIV, 0.9*sigma, 2.0*epsilon, 0.95*reduction); 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); LangevinIntegrator integrator2(0.0, 0.1, 0.01);
Context context2(system, integrator2, Platform::getPlatformByName(platformName)); Context context2(system, integrator2, Platform::getPlatformByName(platformName));
...@@ -195,17 +197,25 @@ void testVdw() { ...@@ -195,17 +197,25 @@ void testVdw() {
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), tolerance); ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), tolerance);
} }
void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff, void setupAndGetForcesEnergyVdwAmmonia2(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff,
double boxDimension, std::vector<Vec3>& forces, double& energy) { double boxDimension, std::vector<Vec3>& forces, double& energy,
AmoebaVdwForce::AlchemicalMethod alchemicalMethod, int softcorePower, double softcoreAlpha, double vdwLambda){
// beginning of Vdw setup // beginning of Vdw setup
System system; System system;
AmoebaVdwForce* amoebaVdwForce = new AmoebaVdwForce();; AmoebaVdwForce* amoebaVdwForce = new AmoebaVdwForce();
int numberOfParticles = 8; int numberOfParticles = 8;
amoebaVdwForce->setSigmaCombiningRule(sigmaCombiningRule); amoebaVdwForce->setSigmaCombiningRule(sigmaCombiningRule);
amoebaVdwForce->setEpsilonCombiningRule(epsilonCombiningRule); amoebaVdwForce->setEpsilonCombiningRule(epsilonCombiningRule);
amoebaVdwForce->setCutoff(cutoff); amoebaVdwForce->setCutoff(cutoff);
bool alchemical = false;
if (alchemicalMethod != AmoebaVdwForce::None) {
amoebaVdwForce->setAlchemicalMethod(alchemicalMethod);
amoebaVdwForce->setSoftcorePower(softcorePower);
amoebaVdwForce->setSoftcoreAlpha(softcoreAlpha);
alchemical = true;
}
if (boxDimension > 0.0) { if (boxDimension > 0.0) {
Vec3 a(boxDimension, 0.0, 0.0); Vec3 a(boxDimension, 0.0, 0.0);
Vec3 b(0.0, boxDimension, 0.0); Vec3 b(0.0, boxDimension, 0.0);
...@@ -213,11 +223,13 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co ...@@ -213,11 +223,13 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co
system.setDefaultPeriodicBoxVectors(a, b, c); system.setDefaultPeriodicBoxVectors(a, b, c);
amoebaVdwForce->setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic); amoebaVdwForce->setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic);
amoebaVdwForce->setUseDispersionCorrection(1); amoebaVdwForce->setUseDispersionCorrection(1);
} else { }
else {
amoebaVdwForce->setNonbondedMethod(AmoebaVdwForce::NoCutoff); amoebaVdwForce->setNonbondedMethod(AmoebaVdwForce::NoCutoff);
amoebaVdwForce->setUseDispersionCorrection(0); amoebaVdwForce->setUseDispersionCorrection(0);
} }
// addParticle: ivIndex, radius, epsilon, reductionFactor // addParticle: ivIndex, radius, epsilon, reductionFactor
system.addParticle( 1.4007000e+01); system.addParticle( 1.4007000e+01);
...@@ -233,16 +245,16 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co ...@@ -233,16 +245,16 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co
amoebaVdwForce->addParticle(0, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01); amoebaVdwForce->addParticle(0, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01);
system.addParticle( 1.4007000e+01); system.addParticle( 1.4007000e+01);
amoebaVdwForce->addParticle(4, 1.8550000e-01, 4.3932000e-01, 0.0000000e+00); amoebaVdwForce->addParticle(4, 1.8550000e-01, 4.3932000e-01, 0.0000000e+00, alchemical);
system.addParticle( 1.0080000e+00); system.addParticle( 1.0080000e+00);
amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01); amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01, alchemical);
system.addParticle( 1.0080000e+00); system.addParticle( 1.0080000e+00);
amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01); amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01, alchemical);
system.addParticle( 1.0080000e+00); system.addParticle( 1.0080000e+00);
amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01); amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01, alchemical);
// ParticleExclusions // ParticleExclusions
...@@ -323,12 +335,28 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co ...@@ -323,12 +335,28 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName(platformName)); Context context(system, integrator, Platform::getPlatformByName(platformName));
// Load the vdw lambda value into the context.
if (alchemicalMethod != AmoebaVdwForce::None) {
context.setParameter(AmoebaVdwForce::Lambda(), vdwLambda);
}
context.setPositions(positions); context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
forces = state.getForces(); forces = state.getForces();
energy = state.getPotentialEnergy(); energy = state.getPotentialEnergy();
} }
void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff,
double boxDimension, std::vector<Vec3>& forces, double& energy) {
AmoebaVdwForce::AlchemicalMethod alchemicalMethod = AmoebaVdwForce::None;
int softcorePower = 5;
double softcoreAlpha = 0.7;
double vdwLambda = 1.0;
setupAndGetForcesEnergyVdwAmmonia2(sigmaCombiningRule, epsilonCombiningRule, cutoff, boxDimension, forces, energy,
alchemicalMethod, softcorePower, softcoreAlpha, vdwLambda);
}
void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy, void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces, std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance) { std::vector<Vec3>& forces, double tolerance) {
...@@ -368,6 +396,43 @@ void testVdwAmmoniaCubicMeanHhg() { ...@@ -368,6 +396,43 @@ void testVdwAmmoniaCubicMeanHhg() {
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
// test alchemical VDW
void testVdwAlchemical(int power, double alpha, double lambda, AmoebaVdwForce::AlchemicalMethod method) {
std::string testName = "testVdwAlchemical";
int numberOfParticles = 8;
double boxDimension = -1.0;
double cutoff = 9000000.0;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyVdwAmmonia2("CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy,
method, power, alpha, lambda);
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 4.8012258e+00;
expectedForces[0] = Vec3( 2.9265247e+02, -1.4507808e-02, -6.9562123e+00);
expectedForces[1] = Vec3(-2.2451693e+00, 4.8143073e-01, -2.0041494e-01);
expectedForces[2] = Vec3(-2.2440698e+00, -4.7905450e-01, -2.0125284e-01);
expectedForces[3] = Vec3(-1.0840394e+00, -5.8531253e-04, 2.6934135e-01);
expectedForces[4] = Vec3(-5.6305662e+01, 1.4733908e-03, -1.8083306e-01);
expectedForces[5] = Vec3( 1.6750145e+00, -3.2448374e-01, -1.8030914e-01);
expectedForces[6] = Vec3(-2.3412420e+02, 1.0754069e-02, 7.6287492e+00);
expectedForces[7] = Vec3( 1.6756544e+00, 3.2497316e-01, -1.7906832e-01);
double scale = pow(lambda, power);
expectedEnergy *= scale;
for (int i=0; i<8; i++) {
expectedForces[i] *= scale;
}
double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
}
// test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic // test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
void testVdwAmmoniaArithmeticArithmetic() { void testVdwAmmoniaArithmeticArithmetic() {
...@@ -2062,6 +2127,25 @@ int main(int argc, char* argv[]) { ...@@ -2062,6 +2127,25 @@ int main(int argc, char* argv[]) {
testTriclinic(); testTriclinic();
// Set lambda and the softcore power (n) to any values, while softcore alpha set to 0.
// The energy and forces are equal to scaling those from testVdwAmmoniaCubicMeanHhg by lambda^n;
int n = 5;
double alpha = 0.0;
double lambda = 0.9;
AmoebaVdwForce::AlchemicalMethod method = AmoebaVdwForce::Annihilate;
testVdwAlchemical(n, alpha, lambda, method);
// Test the Decouple alchemical method.
lambda = 0.5;
method = AmoebaVdwForce::Decouple;
testVdwAlchemical(n, alpha, lambda, method);
// Test alpha > 0.
// This requires lambda = 0, since the energy and forces are not simply scaled by lambda^n.
lambda = 0.0;
alpha = 0.7;
testVdwAlchemical(n, alpha, lambda, method);
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl; std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ 68424035
...@@ -1017,14 +1017,16 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A ...@@ -1017,14 +1017,16 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
sigmas.resize(numParticles); sigmas.resize(numParticles);
epsilons.resize(numParticles); epsilons.resize(numParticles);
reductions.resize(numParticles); reductions.resize(numParticles);
isAlchemical.resize(numParticles);
for (int ii = 0; ii < numParticles; ii++) { for (int ii = 0; ii < numParticles; ii++) {
int indexIV; int indexIV;
double sigma, epsilon, reduction; double sigma, epsilon, reduction;
bool alchemical;
std::vector<int> exclusions; std::vector<int> exclusions;
force.getParticleParameters(ii, indexIV, sigma, epsilon, reduction); force.getParticleParameters(ii, indexIV, sigma, epsilon, reduction, alchemical);
force.getParticleExclusions(ii, exclusions); force.getParticleExclusions(ii, exclusions);
for (unsigned int jj = 0; jj < exclusions.size(); jj++) { for (unsigned int jj = 0; jj < exclusions.size(); jj++) {
allExclusions[ii].insert(exclusions[jj]); allExclusions[ii].insert(exclusions[jj]);
...@@ -1034,6 +1036,7 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A ...@@ -1034,6 +1036,7 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
sigmas[ii] = sigma; sigmas[ii] = sigma;
epsilons[ii] = epsilon; epsilons[ii] = epsilon;
reductions[ii] = reduction; reductions[ii] = reduction;
isAlchemical[ii] = alchemical;
} }
sigmaCombiningRule = force.getSigmaCombiningRule(); sigmaCombiningRule = force.getSigmaCombiningRule();
epsilonCombiningRule = force.getEpsilonCombiningRule(); epsilonCombiningRule = force.getEpsilonCombiningRule();
...@@ -1042,7 +1045,9 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A ...@@ -1042,7 +1045,9 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
cutoff = force.getCutoffDistance(); cutoff = force.getCutoffDistance();
neighborList = useCutoff ? new NeighborList() : NULL; neighborList = useCutoff ? new NeighborList() : NULL;
dispersionCoefficient = force.getUseDispersionCorrection() ? AmoebaVdwForceImpl::calcDispersionCorrection(system, force) : 0.0; dispersionCoefficient = force.getUseDispersionCorrection() ? AmoebaVdwForceImpl::calcDispersionCorrection(system, force) : 0.0;
alchemicalMethod = force.getAlchemicalMethod();
n = force.getSoftcorePower();
alpha = force.getSoftcoreAlpha();
} }
double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) { double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
...@@ -1050,7 +1055,17 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc ...@@ -1050,7 +1055,17 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
vector<Vec3>& posData = extractPositions(context); vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context); vector<Vec3>& forceData = extractForces(context);
AmoebaReferenceVdwForce vdwForce(sigmaCombiningRule, epsilonCombiningRule); AmoebaReferenceVdwForce vdwForce(sigmaCombiningRule, epsilonCombiningRule);
if (alchemicalMethod == AmoebaVdwForce::Decouple) {
vdwForce.setAlchemicalMethod(AmoebaReferenceVdwForce::Decouple);
} else if (alchemicalMethod == AmoebaVdwForce::Annihilate) {
vdwForce.setAlchemicalMethod(AmoebaReferenceVdwForce::Annihilate);
} else {
vdwForce.setAlchemicalMethod(AmoebaReferenceVdwForce::None);
}
vdwForce.setSoftcorePower(n);
vdwForce.setSoftcoreAlpha(alpha);
double energy; double energy;
double lambda = context.getParameter(AmoebaVdwForce::Lambda());
if (useCutoff) { if (useCutoff) {
vdwForce.setCutoff(cutoff); vdwForce.setCutoff(cutoff);
computeNeighborListVoxelHash(*neighborList, numParticles, posData, allExclusions, extractBoxVectors(context), usePBC, cutoff, 0.0); computeNeighborListVoxelHash(*neighborList, numParticles, posData, allExclusions, extractBoxVectors(context), usePBC, cutoff, 0.0);
...@@ -1062,14 +1077,16 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc ...@@ -1062,14 +1077,16 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
throw OpenMMException("The periodic box size has decreased to less than twice the cutoff."); throw OpenMMException("The periodic box size has decreased to less than twice the cutoff.");
} }
vdwForce.setPeriodicBox(boxVectors); 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]); energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
} else { } else {
vdwForce.setNonbondedMethod(AmoebaReferenceVdwForce::CutoffNonPeriodic); vdwForce.setNonbondedMethod(AmoebaReferenceVdwForce::CutoffNonPeriodic);
} }
} else { } else {
vdwForce.setNonbondedMethod(AmoebaReferenceVdwForce::NoCutoff); 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); return static_cast<double>(energy);
} }
...@@ -1079,15 +1096,16 @@ void ReferenceCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& con ...@@ -1079,15 +1096,16 @@ void ReferenceCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& con
throw OpenMMException("updateParametersInContext: The number of particles has changed"); throw OpenMMException("updateParametersInContext: The number of particles has changed");
// Record the values. // Record the values.
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
int indexIV; int indexIV;
double sigma, epsilon, reduction; double sigma, epsilon, reduction;
force.getParticleParameters(i, indexIV, sigma, epsilon, reduction); bool alchemical;
force.getParticleParameters(i, indexIV, sigma, epsilon, reduction, alchemical);
indexIVs[i] = indexIV; indexIVs[i] = indexIV;
sigmas[i] = sigma; sigmas[i] = sigma;
epsilons[i] = epsilon; epsilons[i] = epsilon;
reductions[i]= reduction; reductions[i]= reduction;
isAlchemical[i]= alchemical;
} }
} }
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ 68424035
...@@ -501,11 +501,15 @@ private: ...@@ -501,11 +501,15 @@ private:
int usePBC; int usePBC;
double cutoff; double cutoff;
double dispersionCoefficient; double dispersionCoefficient;
AmoebaVdwForce::AlchemicalMethod alchemicalMethod;
int n;
double alpha;
std::vector<int> indexIVs; std::vector<int> indexIVs;
std::vector< std::set<int> > allExclusions; std::vector< std::set<int> > allExclusions;
std::vector<double> sigmas; std::vector<double> sigmas;
std::vector<double> epsilons; std::vector<double> epsilons;
std::vector<double> reductions; std::vector<double> reductions;
std::vector<bool> isAlchemical;
std::string sigmaCombiningRule; std::string sigmaCombiningRule;
std::string epsilonCombiningRule; std::string epsilonCombiningRule;
const System& system; const System& system;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp
View file @ 68424035
...@@ -32,7 +32,7 @@ ...@@ -32,7 +32,7 @@
using std::vector; using std::vector;
using namespace OpenMM; 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); setTaperCoefficients(_cutoff);
setSigmaCombiningRule("ARITHMETIC"); setSigmaCombiningRule("ARITHMETIC");
...@@ -40,7 +40,7 @@ AmoebaReferenceVdwForce::AmoebaReferenceVdwForce() : _nonbondedMethod(NoCutoff), ...@@ -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); setTaperCoefficients(_cutoff);
setSigmaCombiningRule(sigmaCombiningRule); setSigmaCombiningRule(sigmaCombiningRule);
...@@ -55,6 +55,32 @@ void AmoebaReferenceVdwForce::setNonbondedMethod(AmoebaReferenceVdwForce::Nonbon ...@@ -55,6 +55,32 @@ void AmoebaReferenceVdwForce::setNonbondedMethod(AmoebaReferenceVdwForce::Nonbon
_nonbondedMethod = nonbondedMethod; _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) { void AmoebaReferenceVdwForce::setTaperCoefficients(double cutoff) {
_taperCutoff = cutoff*_taperCutoffFactor; _taperCutoff = cutoff*_taperCutoffFactor;
if (_taperCutoff != cutoff) { if (_taperCutoff != cutoff) {
...@@ -170,13 +196,15 @@ void AmoebaReferenceVdwForce::addReducedForce(unsigned int particleI, unsigned i ...@@ -170,13 +196,15 @@ void AmoebaReferenceVdwForce::addReducedForce(unsigned int particleI, unsigned i
forces[particleIV][2] += sign*force[2]*(1.0 - reduction); 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& particleIPosition,
const Vec3& particleJPosition, const Vec3& particleJPosition,
Vec3& force) const { Vec3& force) const {
static const double dhal = 0.07; static const double dhal = 0.07;
static const double ghal = 0.12; 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 // get deltaR, R2, and R between 2 atoms
...@@ -188,29 +216,25 @@ double AmoebaReferenceVdwForce::calculatePairIxn(double combinedSigma, double co ...@@ -188,29 +216,25 @@ double AmoebaReferenceVdwForce::calculatePairIxn(double combinedSigma, double co
double r_ij_2 = deltaR[ReferenceForce::R2Index]; double r_ij_2 = deltaR[ReferenceForce::R2Index];
double r_ij = deltaR[ReferenceForce::RIndex]; 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 rho = r_ij / combinedSigma;
double r_ij_7 = r_ij_6*r_ij; double rho2 = rho * rho;
double rho6 = rho2 * rho2 * rho2;
double rho = r_ij_7 + ghal*sigma_7; double rhoplus = rho + dhal;
double rhodec2 = rhoplus * rhoplus;
double tau = (dhal + 1.0)/(r_ij + dhal*combinedSigma); double rhodec = rhodec2 * rhodec2 * rhodec2;
double tau_7 = tau*tau*tau; double s1 = 1.0 / (softcore + rhodec * rhoplus);
tau_7 = tau_7*tau_7*tau; double s2 = 1.0 / (softcore + rho6 * rho + 0.12);
double point72 = dhal1 * dhal1;
double dtau = tau/(dhal + 1.0); double t1 = dhal1 * point72 * point72 * point72 * s1;
double t2 = ghal1 * s2;
double ratio = (sigma_7/rho); double t2min = t2 - 2;
double gtau = combinedEpsilon*tau_7*r_ij_6*(ghal+1.0)*ratio*ratio; double dt1 = -7.0 * rhodec * t1 * s1;
double dt2 = -7.0 * rho6 * t2 * s2;
double energy = combinedEpsilon*tau_7*sigma_7*((ghal+1.0)*sigma_7/rho - 2.0); double energy = combinedEpsilon * t1 * t2min;
double dEdR = combinedEpsilon * (dt1 * t2min + t1 * dt2) / combinedSigma;
double dEdR = -7.0*(dtau*energy + gtau);
// tapering // tapering
if ((_nonbondedMethod == CutoffNonPeriodic || _nonbondedMethod == CutoffPeriodic) && r_ij > _taperCutoff) { if ((_nonbondedMethod == CutoffNonPeriodic || _nonbondedMethod == CutoffPeriodic) && r_ij > _taperCutoff) {
double delta = r_ij - _taperCutoff; double delta = r_ij - _taperCutoff;
double taper = 1.0 + delta*delta*delta*(_taperCoefficients[C3] + delta*(_taperCoefficients[C4] + delta*_taperCoefficients[C5])); double taper = 1.0 + delta*delta*delta*(_taperCoefficients[C3] + delta*(_taperCoefficients[C4] + delta*_taperCoefficients[C5]));
...@@ -248,12 +272,13 @@ void AmoebaReferenceVdwForce::setReducedPositions(int numParticles, ...@@ -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 vector<Vec3>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<int>& indexIVs,
const std::vector<double>& sigmas, const std::vector<double>& sigmas,
const std::vector<double>& epsilons, const std::vector<double>& epsilons,
const std::vector<double>& reductions, const std::vector<double>& reductions,
const std::vector<bool>& isAlchemical,
const std::vector< std::set<int> >& allExclusions, const std::vector< std::set<int> >& allExclusions,
vector<Vec3>& forces) const { vector<Vec3>& forces) const {
...@@ -277,6 +302,7 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, ...@@ -277,6 +302,7 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
double sigmaI = sigmas[ii]; double sigmaI = sigmas[ii];
double epsilonI = epsilons[ii]; double epsilonI = epsilons[ii];
bool isAlchemicalI = isAlchemical[ii];
for (int jj : allExclusions[ii]) for (int jj : allExclusions[ii])
exclusions[jj] = 1; exclusions[jj] = 1;
...@@ -285,12 +311,20 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, ...@@ -285,12 +311,20 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
double combinedSigma = (this->*_combineSigmas)(sigmaI, sigmas[jj]); double combinedSigma = (this->*_combineSigmas)(sigmaI, sigmas[jj]);
double combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj]); double combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj]);
double softcore = 0.0;
if (this->_alchemicalMethod == Decouple && (isAlchemicalI != isAlchemical[jj])) {
combinedEpsilon *= pow(lambda, this->_n);
softcore = this->_alpha * pow(1.0 - lambda, 2);
} else if (this->_alchemicalMethod == Annihilate && (isAlchemicalI || isAlchemical[jj])) {
combinedEpsilon *= pow(lambda, this->_n);
softcore = this->_alpha * pow(1.0 - lambda, 2);
}
Vec3 force; Vec3 force;
energy += calculatePairIxn(combinedSigma, combinedEpsilon, energy += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
reducedPositions[ii], reducedPositions[jj], reducedPositions[ii], reducedPositions[jj], force);
force);
if (indexIVs[ii] == ii) { if (indexIVs[ii] == ii) {
forces[ii][0] -= force[0]; forces[ii][0] -= force[0];
forces[ii][1] -= force[1]; forces[ii][1] -= force[1];
...@@ -316,12 +350,13 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, ...@@ -316,12 +350,13 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
return energy; return energy;
} }
double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double lambda,
const vector<Vec3>& particlePositions, const vector<Vec3>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<int>& indexIVs,
const std::vector<double>& sigmas, const std::vector<double>& sigmas,
const std::vector<double>& epsilons, const std::vector<double>& epsilons,
const std::vector<double>& reductions, const std::vector<double>& reductions,
const std::vector<bool>& isAlchemical,
const NeighborList& neighborList, const NeighborList& neighborList,
vector<Vec3>& forces) const { vector<Vec3>& forces) const {
...@@ -329,7 +364,7 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, ...@@ -329,7 +364,7 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
std::vector<Vec3> reducedPositions; std::vector<Vec3> reducedPositions;
setReducedPositions(numParticles, particlePositions, indexIVs, reductions, reducedPositions); setReducedPositions(numParticles, particlePositions, indexIVs, reductions, reducedPositions);
// loop over neighbor list // loop over neighbor list
// (1) calculate pair vdw ixn // (1) calculate pair vdw ixn
// (2) accumulate forces: if particle is a site where interaction position != particle position, // (2) accumulate forces: if particle is a site where interaction position != particle position,
...@@ -345,10 +380,21 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, ...@@ -345,10 +380,21 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
double combinedSigma = (this->*_combineSigmas)(sigmas[siteI], sigmas[siteJ]); double combinedSigma = (this->*_combineSigmas)(sigmas[siteI], sigmas[siteJ]);
double combinedEpsilon = (this->*_combineEpsilons)(epsilons[siteI], epsilons[siteJ]); double combinedEpsilon = (this->*_combineEpsilons)(epsilons[siteI], epsilons[siteJ]);
double softcore = 0.0;
int isAlchemicalI = isAlchemical[siteI];
int isAlchemicalJ = isAlchemical[siteJ];
if (this->_alchemicalMethod == Decouple && (isAlchemicalI != isAlchemicalJ)) {
combinedEpsilon *= pow(lambda, this->_n);
softcore = this->_alpha * pow(1.0 - lambda, 2);
} else if (this->_alchemicalMethod == Annihilate && (isAlchemicalI || isAlchemicalJ)) {
combinedEpsilon *= pow(lambda, this->_n);
softcore = this->_alpha * pow(1.0 - lambda, 2);
}
Vec3 force; Vec3 force;
energy += calculatePairIxn(combinedSigma, combinedEpsilon, energy += calculatePairIxn(combinedSigma, combinedEpsilon, softcore,
reducedPositions[siteI], reducedPositions[siteJ], force); reducedPositions[siteI], reducedPositions[siteJ], force);
if (indexIVs[siteI] == siteI) { if (indexIVs[siteI] == siteI) {
forces[siteI][0] -= force[0]; forces[siteI][0] -= force[0];
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.h
View file @ 68424035
...@@ -63,6 +63,25 @@ public: ...@@ -63,6 +63,25 @@ public:
*/ */
CutoffPeriodic = 2, 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: ...@@ -109,6 +128,27 @@ public:
void setNonbondedMethod(NonbondedMethod nonbondedMethod); 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 Get cutoff
...@@ -129,6 +169,47 @@ public: ...@@ -129,6 +169,47 @@ public:
void setCutoff(double cutoff); 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 Set sigma combining rule
...@@ -184,11 +265,13 @@ public: ...@@ -184,11 +265,13 @@ public:
Calculate Amoeba Hal vdw ixns Calculate Amoeba Hal vdw ixns
@param numParticles number of particles @param numParticles number of particles
@param lambda lambda value
@param particlePositions Cartesian coordinates of particles @param particlePositions Cartesian coordinates of particles
@param indexIVs position index for associated reducing particle @param indexIVs position index for associated reducing particle
@param sigmas particle sigmas @param sigmas particle sigmas
@param epsilons particle epsilons @param epsilons particle epsilons
@param reductions particle reduction factors @param reductions particle reduction factors
@param isAlchemical particle alchemical flag
@param vdwExclusions particle exclusions @param vdwExclusions particle exclusions
@param forces add forces to this vector @param forces add forces to this vector
...@@ -196,10 +279,11 @@ public: ...@@ -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<int>& indexIVs,
const std::vector<double>& sigmas, const std::vector<double>& epsilons, const std::vector<double>& sigmas, const std::vector<double>& epsilons,
const std::vector<double>& reductions, const std::vector<double>& reductions,
const std::vector<bool>& isAlchemical,
const std::vector< std::set<int> >& vdwExclusions, const std::vector< std::set<int> >& vdwExclusions,
std::vector<OpenMM::Vec3>& forces) const; std::vector<OpenMM::Vec3>& forces) const;
...@@ -208,11 +292,13 @@ public: ...@@ -208,11 +292,13 @@ public:
Calculate Vdw ixn using neighbor list Calculate Vdw ixn using neighbor list
@param numParticles number of particles @param numParticles number of particles
@param lambda lambda value
@param particlePositions Cartesian coordinates of particles @param particlePositions Cartesian coordinates of particles
@param indexIVs position index for associated reducing particle @param indexIVs position index for associated reducing particle
@param sigmas particle sigmas @param sigmas particle sigmas
@param epsilons particle epsilons @param epsilons particle epsilons
@param reductions particle reduction factors @param reductions particle reduction factors
@param isAlchemical particle alchemical flag
@param neighborList neighbor list @param neighborList neighbor list
@param forces add forces to this vector @param forces add forces to this vector
...@@ -220,17 +306,16 @@ public: ...@@ -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<int>& indexIVs,
const std::vector<double>& sigmas, const std::vector<double>& epsilons, const std::vector<double>& sigmas, const std::vector<double>& epsilons,
const std::vector<double>& reductions, const std::vector<double>& reductions,
const std::vector<bool>& isAlchemical,
const NeighborList& neighborList, const NeighborList& neighborList,
std::vector<OpenMM::Vec3>& forces) const; std::vector<OpenMM::Vec3>& forces) const;
private: private:
// taper coefficient indices // taper coefficient indices
static const int C3=0; static const int C3=0;
static const int C4=1; static const int C4=1;
static const int C5=2; static const int C5=2;
...@@ -238,6 +323,9 @@ private: ...@@ -238,6 +323,9 @@ private:
std::string _sigmaCombiningRule; std::string _sigmaCombiningRule;
std::string _epsilonCombiningRule; std::string _epsilonCombiningRule;
NonbondedMethod _nonbondedMethod; NonbondedMethod _nonbondedMethod;
AlchemicalMethod _alchemicalMethod;
double _n;
double _alpha;
double _cutoff; double _cutoff;
double _taperCutoffFactor; double _taperCutoffFactor;
double _taperCutoff; double _taperCutoff;
...@@ -245,14 +333,13 @@ private: ...@@ -245,14 +333,13 @@ private:
Vec3 _periodicBoxVectors[3]; Vec3 _periodicBoxVectors[3];
CombiningFunction _combineSigmas; CombiningFunction _combineSigmas;
double arithmeticSigmaCombiningRule(double sigmaI, double sigmaJ) const; double arithmeticSigmaCombiningRule(double sigmaI, double sigmaJ) const;
double geometricSigmaCombiningRule(double sigmaI, double sigmaJ) const; double geometricSigmaCombiningRule(double sigmaI, double sigmaJ) const;
double cubicMeanSigmaCombiningRule(double sigmaI, double sigmaJ) const; double cubicMeanSigmaCombiningRule(double sigmaI, double sigmaJ) const;
CombiningFunction _combineEpsilons; CombiningFunction _combineEpsilons;
double arithmeticEpsilonCombiningRule(double epsilonI, double epsilonJ) const; double arithmeticEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
double geometricEpsilonCombiningRule(double epsilonI, double epsilonJ) const; double geometricEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
double harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ) const; double harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
double hhgEpsilonCombiningRule(double epsilonI, double epsilonJ) const; double hhgEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -308,6 +395,7 @@ private: ...@@ -308,6 +395,7 @@ private:
@param combindedSigma combined sigmas @param combindedSigma combined sigmas
@param combindedEpsilon combined epsilons @param combindedEpsilon combined epsilons
@param softcore softcore offset parameter
@param particleIPosition particle I position @param particleIPosition particle I position
@param particleJPosition particle J position @param particleJPosition particle J position
@param force output force @param force output force
...@@ -316,7 +404,7 @@ private: ...@@ -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, const Vec3& particleIPosition, const Vec3& particleJPosition,
Vec3& force) const; Vec3& force) const;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaVdwForce.cpp
View file @ 68424035
...@@ -142,10 +142,11 @@ void testVdw() { ...@@ -142,10 +142,11 @@ void testVdw() {
for (int ii = 0; ii < amoebaVdwForce->getNumParticles(); ii++) { for (int ii = 0; ii < amoebaVdwForce->getNumParticles(); ii++) {
int indexIV; int indexIV;
double sigma, epsilon, reduction; double sigma, epsilon, reduction;
amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction); bool isAlchemical;
amoebaVdwForce->getParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
sigma *= AngstromToNm; sigma *= AngstromToNm;
epsilon *= CalToJoule; epsilon *= CalToJoule;
amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction); amoebaVdwForce->setParticleParameters(ii, indexIV, sigma, epsilon, reduction, isAlchemical);
} }
platformName = "Reference"; platformName = "Reference";
Context context(system, integrator, Platform::getPlatformByName(platformName)); Context context(system, integrator, Platform::getPlatformByName(platformName));
...@@ -173,8 +174,9 @@ void testVdw() { ...@@ -173,8 +174,9 @@ void testVdw() {
for (int i = 0; i < numberOfParticles; i++) { for (int i = 0; i < numberOfParticles; i++) {
int indexIV; int indexIV;
double mass, sigma, epsilon, reduction; double mass, sigma, epsilon, reduction;
amoebaVdwForce->getParticleParameters(i, indexIV, sigma, epsilon, reduction); bool isAlchemical;
amoebaVdwForce->setParticleParameters(i, indexIV, 0.9*sigma, 2.0*epsilon, 0.95*reduction); 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); LangevinIntegrator integrator2(0.0, 0.1, 0.01);
Context context2(system, integrator2, Platform::getPlatformByName(platformName)); Context context2(system, integrator2, Platform::getPlatformByName(platformName));
...@@ -198,17 +200,26 @@ void testVdw() { ...@@ -198,17 +200,26 @@ void testVdw() {
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), tolerance); ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), tolerance);
} }
void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff,
double boxDimension, std::vector<Vec3>& forces, double& energy) { void setupAndGetForcesEnergyVdwAmmonia2(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff,
double boxDimension, std::vector<Vec3>& forces, double& energy,
AmoebaVdwForce::AlchemicalMethod alchemicalMethod, int softcorePower, double softcoreAlpha, double vdwLambda){
// beginning of Vdw setup // beginning of Vdw setup
System system; System system;
AmoebaVdwForce* amoebaVdwForce = new AmoebaVdwForce();; AmoebaVdwForce* amoebaVdwForce = new AmoebaVdwForce();
int numberOfParticles = 8; int numberOfParticles = 8;
amoebaVdwForce->setSigmaCombiningRule(sigmaCombiningRule); amoebaVdwForce->setSigmaCombiningRule(sigmaCombiningRule);
amoebaVdwForce->setEpsilonCombiningRule(epsilonCombiningRule); amoebaVdwForce->setEpsilonCombiningRule(epsilonCombiningRule);
amoebaVdwForce->setCutoff(cutoff); amoebaVdwForce->setCutoff(cutoff);
bool alchemical = false;
if (alchemicalMethod != AmoebaVdwForce::None) {
amoebaVdwForce->setAlchemicalMethod(alchemicalMethod);
amoebaVdwForce->setSoftcorePower(softcorePower);
amoebaVdwForce->setSoftcoreAlpha(softcoreAlpha);
alchemical = true;
}
if (boxDimension > 0.0) { if (boxDimension > 0.0) {
Vec3 a(boxDimension, 0.0, 0.0); Vec3 a(boxDimension, 0.0, 0.0);
Vec3 b(0.0, boxDimension, 0.0); Vec3 b(0.0, boxDimension, 0.0);
...@@ -237,16 +248,16 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co ...@@ -237,16 +248,16 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co
amoebaVdwForce->addParticle(0, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01); amoebaVdwForce->addParticle(0, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01);
system.addParticle( 1.4007000e+01); system.addParticle( 1.4007000e+01);
amoebaVdwForce->addParticle(4, 1.8550000e-01, 4.3932000e-01, 0.0000000e+00); amoebaVdwForce->addParticle(4, 1.8550000e-01, 4.3932000e-01, 0.0000000e+00, alchemical);
system.addParticle( 1.0080000e+00); system.addParticle( 1.0080000e+00);
amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01); amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01, alchemical);
system.addParticle( 1.0080000e+00); system.addParticle( 1.0080000e+00);
amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01); amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01, alchemical);
system.addParticle( 1.0080000e+00); system.addParticle( 1.0080000e+00);
amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01); amoebaVdwForce->addParticle(4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01, alchemical);
// ParticleExclusions // ParticleExclusions
...@@ -334,12 +345,27 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co ...@@ -334,12 +345,27 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName(platformName)); Context context(system, integrator, Platform::getPlatformByName(platformName));
// Load the vdw lambda value into the context.
if (alchemicalMethod != AmoebaVdwForce::None) {
context.setParameter(AmoebaVdwForce::Lambda(), vdwLambda);
}
context.setPositions(positions); context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
forces = state.getForces(); forces = state.getForces();
energy = state.getPotentialEnergy(); energy = state.getPotentialEnergy();
} }
void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff,
double boxDimension, std::vector<Vec3>& forces, double& energy) {
AmoebaVdwForce::AlchemicalMethod alchemicalMethod = AmoebaVdwForce::None;
int softcorePower = 5;
double softcoreAlpha = 0.7;
double vdwLambda = 1.0;
setupAndGetForcesEnergyVdwAmmonia2(sigmaCombiningRule, epsilonCombiningRule, cutoff, boxDimension, forces, energy,
alchemicalMethod, softcorePower, softcoreAlpha, vdwLambda);
}
void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy, void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces, std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance) { std::vector<Vec3>& forces, double tolerance) {
...@@ -379,6 +405,43 @@ void testVdwAmmoniaCubicMeanHhg() { ...@@ -379,6 +405,43 @@ void testVdwAmmoniaCubicMeanHhg() {
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
// test alchemical VDW
void testVdwAlchemical(int power, double alpha, double lambda, AmoebaVdwForce::AlchemicalMethod method) {
std::string testName = "testVdwAlchemical";
int numberOfParticles = 8;
double boxDimension = -1.0;
double cutoff = 9000000.0;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyVdwAmmonia2("CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy,
method, power, alpha, lambda);
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 4.8012258e+00;
expectedForces[0] = Vec3( 2.9265247e+02, -1.4507808e-02, -6.9562123e+00);
expectedForces[1] = Vec3(-2.2451693e+00, 4.8143073e-01, -2.0041494e-01);
expectedForces[2] = Vec3(-2.2440698e+00, -4.7905450e-01, -2.0125284e-01);
expectedForces[3] = Vec3(-1.0840394e+00, -5.8531253e-04, 2.6934135e-01);
expectedForces[4] = Vec3(-5.6305662e+01, 1.4733908e-03, -1.8083306e-01);
expectedForces[5] = Vec3( 1.6750145e+00, -3.2448374e-01, -1.8030914e-01);
expectedForces[6] = Vec3(-2.3412420e+02, 1.0754069e-02, 7.6287492e+00);
expectedForces[7] = Vec3( 1.6756544e+00, 3.2497316e-01, -1.7906832e-01);
double scale = pow(lambda, power);
expectedEnergy *= scale;
for (int i=0; i<8; i++) {
expectedForces[i] *= scale;
}
double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
}
// test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic // test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
void testVdwAmmoniaArithmeticArithmetic() { void testVdwAmmoniaArithmeticArithmetic() {
...@@ -2074,6 +2137,25 @@ int main(int numberOfArguments, char* argv[]) { ...@@ -2074,6 +2137,25 @@ int main(int numberOfArguments, char* argv[]) {
testTriclinic(); testTriclinic();
// Set lambda and the softcore power (n) to any values (softcore alpha set to 0).
// The energy and forces are equal to scaling testVdwAmmoniaCubicMeanHhg by lambda^n;
int n = 5;
double alpha = 0.0;
double lambda = 0.9;
AmoebaVdwForce::AlchemicalMethod method = AmoebaVdwForce::Annihilate;
testVdwAlchemical(n, alpha, lambda, method);
// Test the Decouple alchemical method.
lambda = 0.5;
method = AmoebaVdwForce::Decouple;
testVdwAlchemical(n, alpha, lambda, method);
// Test alpha > 0.
// This requires lambda = 0, since the energy and forces are not simply scaled by lambda^n.
lambda = 0.0;
alpha = 0.7;
testVdwAlchemical(n, alpha, lambda, method);
} }
catch(const std::exception& e) { catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/serialization/src/AmoebaVdwForceProxy.cpp
View file @ 68424035
...@@ -42,25 +42,29 @@ AmoebaVdwForceProxy::AmoebaVdwForceProxy() : SerializationProxy("AmoebaVdwForce" ...@@ -42,25 +42,29 @@ AmoebaVdwForceProxy::AmoebaVdwForceProxy() : SerializationProxy("AmoebaVdwForce"
} }
void AmoebaVdwForceProxy::serialize(const void* object, SerializationNode& node) const { void AmoebaVdwForceProxy::serialize(const void* object, SerializationNode& node) const {
node.setIntProperty("version", 2); node.setIntProperty("version", 3);
const AmoebaVdwForce& force = *reinterpret_cast<const AmoebaVdwForce*>(object); const AmoebaVdwForce& force = *reinterpret_cast<const AmoebaVdwForce*>(object);
node.setIntProperty("forceGroup", force.getForceGroup()); node.setIntProperty("forceGroup", force.getForceGroup());
node.setStringProperty("SigmaCombiningRule", force.getSigmaCombiningRule()); node.setStringProperty("SigmaCombiningRule", force.getSigmaCombiningRule());
node.setStringProperty("EpsilonCombiningRule", force.getEpsilonCombiningRule()); node.setStringProperty("EpsilonCombiningRule", force.getEpsilonCombiningRule());
node.setDoubleProperty("VdwCutoff", force.getCutoffDistance()); node.setDoubleProperty("VdwCutoff", force.getCutoffDistance());
node.setIntProperty("method", (int) force.getNonbondedMethod()); 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"); SerializationNode& particles = node.createChildNode("VdwParticles");
for (unsigned int ii = 0; ii < static_cast<unsigned int>(force.getNumParticles()); ii++) { for (unsigned int ii = 0; ii < static_cast<unsigned int>(force.getNumParticles()); ii++) {
int ivIndex; int ivIndex;
double sigma, epsilon, reductionFactor; 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"); 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; std::vector< int > exclusions;
force.getParticleExclusions(ii, exclusions); force.getParticleExclusions(ii, exclusions);
...@@ -74,7 +78,7 @@ void AmoebaVdwForceProxy::serialize(const void* object, SerializationNode& node) ...@@ -74,7 +78,7 @@ void AmoebaVdwForceProxy::serialize(const void* object, SerializationNode& node)
void* AmoebaVdwForceProxy::deserialize(const SerializationNode& node) const { void* AmoebaVdwForceProxy::deserialize(const SerializationNode& node) const {
int version = node.getIntProperty("version"); int version = node.getIntProperty("version");
if (version < 1 || version > 2) if (version < 1 || version > 3)
throw OpenMMException("Unsupported version number"); throw OpenMMException("Unsupported version number");
AmoebaVdwForce* force = new AmoebaVdwForce(); AmoebaVdwForce* force = new AmoebaVdwForce();
try { try {
...@@ -85,10 +89,23 @@ void* AmoebaVdwForceProxy::deserialize(const SerializationNode& node) const { ...@@ -85,10 +89,23 @@ void* AmoebaVdwForceProxy::deserialize(const SerializationNode& node) const {
force->setCutoffDistance(node.getDoubleProperty("VdwCutoff")); force->setCutoffDistance(node.getDoubleProperty("VdwCutoff"));
force->setNonbondedMethod((AmoebaVdwForce::NonbondedMethod) node.getIntProperty("method")); force->setNonbondedMethod((AmoebaVdwForce::NonbondedMethod) node.getIntProperty("method"));
if (version > 2) {
force->setAlchemicalMethod((AmoebaVdwForce::AlchemicalMethod) node.getIntProperty("alchemicalMethod"));
force->setSoftcorePower(node.getDoubleProperty("n"));
force->setSoftcoreAlpha(node.getDoubleProperty("alpha"));
}
const SerializationNode& particles = node.getChildNode("VdwParticles"); const SerializationNode& particles = node.getChildNode("VdwParticles");
for (unsigned int ii = 0; ii < particles.getChildren().size(); ii++) { for (unsigned int ii = 0; ii < particles.getChildren().size(); ii++) {
const SerializationNode& particle = particles.getChildren()[ii]; const SerializationNode& particle = particles.getChildren()[ii];
force->addParticle(particle.getIntProperty("ivIndex"), particle.getDoubleProperty("sigma"), particle.getDoubleProperty("epsilon"), particle.getDoubleProperty("reductionFactor"));
if (version < 3)
force->addParticle(particle.getIntProperty("ivIndex"), particle.getDoubleProperty("sigma"),
particle.getDoubleProperty("epsilon"), particle.getDoubleProperty("reductionFactor"));
else
force->addParticle(particle.getIntProperty("ivIndex"), particle.getDoubleProperty("sigma"),
particle.getDoubleProperty("epsilon"), particle.getDoubleProperty("reductionFactor"),
particle.getBoolProperty("isAlchemical"));
// exclusions // exclusions
......
plugins/amoeba/serialization/tests/TestSerializeAmoebaVdwForce.cpp
View file @ 68424035
...@@ -50,10 +50,11 @@ void testSerialization() { ...@@ -50,10 +50,11 @@ void testSerialization() {
force1.setEpsilonCombiningRule("GEOMETRIC"); force1.setEpsilonCombiningRule("GEOMETRIC");
force1.setCutoff(0.9); force1.setCutoff(0.9);
force1.setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic); force1.setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic);
force1.setAlchemicalMethod(AmoebaVdwForce::None);
force1.addParticle(0, 1.0, 2.0, 0.9); force1.addParticle(0, 1.0, 2.0, 0.9, false);
force1.addParticle(1, 1.1, 2.1, 0.9); force1.addParticle(1, 1.1, 2.1, 0.9, true);
force1.addParticle(2, 1.3, 4.1, 0.9); force1.addParticle(2, 1.3, 4.1, 0.9, false);
for (unsigned int ii = 0; ii < 3; ii++) { for (unsigned int ii = 0; ii < 3; ii++) {
std::vector< int > exclusions; std::vector< int > exclusions;
exclusions.push_back(ii); exclusions.push_back(ii);
...@@ -76,6 +77,7 @@ void testSerialization() { ...@@ -76,6 +77,7 @@ void testSerialization() {
ASSERT_EQUAL(force1.getEpsilonCombiningRule(), force2.getEpsilonCombiningRule()); ASSERT_EQUAL(force1.getEpsilonCombiningRule(), force2.getEpsilonCombiningRule());
ASSERT_EQUAL(force1.getCutoff(), force2.getCutoff()); ASSERT_EQUAL(force1.getCutoff(), force2.getCutoff());
ASSERT_EQUAL(force1.getNonbondedMethod(), force2.getNonbondedMethod()); ASSERT_EQUAL(force1.getNonbondedMethod(), force2.getNonbondedMethod());
ASSERT_EQUAL(force1.getAlchemicalMethod(), force2.getAlchemicalMethod());
ASSERT_EQUAL(force1.getNumParticles(), force2.getNumParticles()); ASSERT_EQUAL(force1.getNumParticles(), force2.getNumParticles());
...@@ -87,13 +89,17 @@ void testSerialization() { ...@@ -87,13 +89,17 @@ void testSerialization() {
double sigma1, epsilon1, reductionFactor1; double sigma1, epsilon1, reductionFactor1;
double sigma2, epsilon2, reductionFactor2; double sigma2, epsilon2, reductionFactor2;
force1.getParticleParameters(ii, ivIndex1, sigma1, epsilon1, reductionFactor1); bool isAlchemical1;
force2.getParticleParameters(ii, ivIndex2, sigma2, epsilon2, reductionFactor2); bool isAlchemical2;
force1.getParticleParameters(ii, ivIndex1, sigma1, epsilon1, reductionFactor1, isAlchemical1);
force2.getParticleParameters(ii, ivIndex2, sigma2, epsilon2, reductionFactor2, isAlchemical2);
ASSERT_EQUAL(ivIndex1, ivIndex2); ASSERT_EQUAL(ivIndex1, ivIndex2);
ASSERT_EQUAL(sigma1, sigma2); ASSERT_EQUAL(sigma1, sigma2);
ASSERT_EQUAL(epsilon1, epsilon2); ASSERT_EQUAL(epsilon1, epsilon2);
ASSERT_EQUAL(reductionFactor1, reductionFactor2); ASSERT_EQUAL(reductionFactor1, reductionFactor2);
ASSERT_EQUAL(isAlchemical1, isAlchemical2);
} }
for (unsigned int ii = 0; ii < static_cast<unsigned int>(force1.getNumParticles()); ii++) { for (unsigned int ii = 0; ii < static_cast<unsigned int>(force1.getNumParticles()); ii++) {
......
wrappers/python/src/swig_doxygen/swigInputConfig.py
View file @ 68424035
...@@ -309,15 +309,14 @@ UNITS = { ...@@ -309,15 +309,14 @@ UNITS = {
("AmoebaTorsionTorsionForce", "getTorsionTorsionParameters") : ( None, ()), ("AmoebaTorsionTorsionForce", "getTorsionTorsionParameters") : ( None, ()),
("AmoebaTorsionTorsionForce", "getTorsionTorsionGrid") : ( None, ()), ("AmoebaTorsionTorsionForce", "getTorsionTorsionGrid") : ( None, ()),
# LPW 2012-10 : Is this a duplicate entry?
#("AmoebaVdwForce", "getParticleParameters") : ( None, (None, None, None, 'unit.nanometer', 'unit.kilojoule_per_mole', None)),
("AmoebaVdwForce", "getSigmaCombiningRule") : ( None, ()), ("AmoebaVdwForce", "getSigmaCombiningRule") : ( None, ()),
("AmoebaVdwForce", "getEpsilonCombiningRule") : ( None, ()), ("AmoebaVdwForce", "getEpsilonCombiningRule") : ( None, ()),
("AmoebaVdwForce", "getParticleExclusions") : ( None, ()), ("AmoebaVdwForce", "getParticleExclusions") : ( None, ()),
("AmoebaVdwForce", "getAlchemicalMethod") : ( None, ()),
("AmoebaVdwForce", "getSoftcorePower") : ( None, ()),
("AmoebaVdwForce", "getSoftcoreAlpha") : ( None, ()),
("AmoebaVdwForce", "getCutoff") : ( 'unit.nanometer', ()), ("AmoebaVdwForce", "getCutoff") : ( 'unit.nanometer', ()),
# LPW 2012-10 Modified because it no longer returns ivIndex and classIndex. ("AmoebaVdwForce", "getParticleParameters") : ( None, (None, 'unit.nanometer', 'unit.kilojoule_per_mole', None, None)),
# ("AmoebaVdwForce", "getParticleParameters") : ( None, (None, None, 'unit.nanometer', 'unit.kilojoule_per_mole', None)),
("AmoebaVdwForce", "getParticleParameters") : ( None, (None, 'unit.nanometer', 'unit.kilojoule_per_mole', None)),
("AmoebaWcaDispersionForce", "getParticleParameters") : ( None, ('unit.nanometer', 'unit.kilojoule_per_mole')), ("AmoebaWcaDispersionForce", "getParticleParameters") : ( None, ('unit.nanometer', 'unit.kilojoule_per_mole')),
("AmoebaWcaDispersionForce", "getAwater") : ( '1/(unit.nanometer*unit.nanometer*unit.nanometer)',()), ("AmoebaWcaDispersionForce", "getAwater") : ( '1/(unit.nanometer*unit.nanometer*unit.nanometer)',()),
......
wrappers/python/tests/TestAPIUnits.py
View file @ 68424035
...@@ -969,7 +969,7 @@ class TestAPIUnits(unittest.TestCase): ...@@ -969,7 +969,7 @@ class TestAPIUnits(unittest.TestCase):
self.assertEqual(force.getNumParticles(), 3) self.assertEqual(force.getNumParticles(), 3)
p, sig, eps, scale = force.getParticleParameters(0) p, sig, eps, scale, alchemical = force.getParticleParameters(0)
self.assertEqual(p, 0) self.assertEqual(p, 0)
self.assertEqual(sig, 0.1*nanometers) self.assertEqual(sig, 0.1*nanometers)
self.assertIs(sig.unit, nanometers) self.assertIs(sig.unit, nanometers)
...@@ -977,7 +977,7 @@ class TestAPIUnits(unittest.TestCase): ...@@ -977,7 +977,7 @@ class TestAPIUnits(unittest.TestCase):
self.assertIs(eps.unit, kilojoules_per_mole) self.assertIs(eps.unit, kilojoules_per_mole)
self.assertEqual(scale, 1.0) self.assertEqual(scale, 1.0)
p, sig, eps, scale = force.getParticleParameters(1) p, sig, eps, scale, alchemical = force.getParticleParameters(1)
self.assertEqual(p, 1) self.assertEqual(p, 1)
self.assertEqual(sig, 1.0*angstroms) self.assertEqual(sig, 1.0*angstroms)
self.assertIs(sig.unit, nanometers) self.assertIs(sig.unit, nanometers)
...@@ -985,7 +985,7 @@ class TestAPIUnits(unittest.TestCase): ...@@ -985,7 +985,7 @@ class TestAPIUnits(unittest.TestCase):
self.assertIs(eps.unit, kilojoules_per_mole) self.assertIs(eps.unit, kilojoules_per_mole)
self.assertEqual(scale, 0.5) self.assertEqual(scale, 0.5)
p, sig, eps, scale = force.getParticleParameters(2) p, sig, eps, scale, alchemical = force.getParticleParameters(2)
self.assertEqual(p, 1) self.assertEqual(p, 1)
self.assertAlmostEqualUnit(sig, 0.8*angstroms) self.assertAlmostEqualUnit(sig, 0.8*angstroms)
self.assertIs(sig.unit, nanometers) self.assertIs(sig.unit, nanometers)
......
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