Skip to content

GitLab

Commit 9c95f25c authored by Michael Schnieders's avatar Michael Schnieders
Browse files

Added Amoeba vdW softcore tests to both the reference and Cuda platforms

parent 0e55d0f8
Hide whitespace changes
Inline Side-by-side
Showing with 196 additions and 29 deletions
plugins/amoeba/platforms/cuda/src/kernels/amoebaVdwForce2.cu
View file @ 9c95f25c
......@@ -29,13 +29,13 @@
real softcore = 0.0f;
#if VDW_ALCHEMICAL_METHOD == 1
if (isAlchemical1 != isAlchemical2) {
epsilon = epsilon * pow(VDW_LAMBDA, VDW_SOFTCORE_POWER);
softcore = VDW_SOFTCORE_ALPHA * (1.0 - vdwLambda) * (1.0 - vdwLambda);
}
#elif VDW_ALCHEMICAL_METHOD == 2
if (isAlchemical1 || isAlchemical2) {
epsilon = epsilon * pow(VDW_LAMBDA, VDW_SOFTCORE_POWER);
softcore = VDW_SOFTCORE_ALPHA * (1.0 - vdwLambda) * (1.0 - vdwLambda);
#endif
#if VDW_ALCHEMICAL_METHOD != 0
real lambda = vdwLambda[0];
epsilon = epsilon * POW(lambda, VDW_SOFTCORE_POWER);
softcore = VDW_SOFTCORE_ALPHA * (1.0f - lambda) * (1.0f - lambda);
}
#endif
real dhal = 0.07f;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaVdwForce.cpp
View file @ 9c95f25c
......@@ -197,17 +197,25 @@ void testVdw() {
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
System system;
AmoebaVdwForce* amoebaVdwForce = new AmoebaVdwForce();;
AmoebaVdwForce* amoebaVdwForce = new AmoebaVdwForce();
int numberOfParticles = 8;
amoebaVdwForce->setSigmaCombiningRule(sigmaCombiningRule);
amoebaVdwForce->setEpsilonCombiningRule(epsilonCombiningRule);
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) {
Vec3 a(boxDimension, 0.0, 0.0);
Vec3 b(0.0, boxDimension, 0.0);
......@@ -215,11 +223,13 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co
system.setDefaultPeriodicBoxVectors(a, b, c);
amoebaVdwForce->setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic);
amoebaVdwForce->setUseDispersionCorrection(1);
} else {
}
else {
amoebaVdwForce->setNonbondedMethod(AmoebaVdwForce::NoCutoff);
amoebaVdwForce->setUseDispersionCorrection(0);
}
// addParticle: ivIndex, radius, epsilon, reductionFactor
system.addParticle( 1.4007000e+01);
......@@ -235,16 +245,16 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co
amoebaVdwForce->addParticle(0, 1.3500000e-01, 8.3680000e-02, 9.1000000e-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);
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);
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);
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
......@@ -325,12 +335,28 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co
LangevinIntegrator integrator(0.0, 0.1, 0.01);
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);
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
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,
std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance) {
......@@ -370,6 +396,43 @@ void testVdwAmmoniaCubicMeanHhg() {
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
void testVdwAmmoniaArithmeticArithmetic() {
......@@ -2064,6 +2127,25 @@ int main(int argc, char* argv[]) {
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) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ 9c95f25c
......@@ -996,9 +996,6 @@ ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(std::string
CalcAmoebaVdwForceKernel(name, platform), system(system) {
useCutoff = 0;
usePBC = 0;
alchemicalMethod = 0;
n = 5;
alpha = 0.7;
cutoff = 1.0e+10;
neighborList = NULL;
}
......@@ -1048,7 +1045,7 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
cutoff = force.getCutoffDistance();
neighborList = useCutoff ? new NeighborList() : NULL;
dispersionCoefficient = force.getUseDispersionCorrection() ? AmoebaVdwForceImpl::calcDispersionCorrection(system, force) : 0.0;
alchemicalMethod = (int) force.getAlchemicalMethod();
alchemicalMethod = force.getAlchemicalMethod();
n = force.getSoftcorePower();
alpha = force.getSoftcoreAlpha();
}
......@@ -1058,6 +1055,15 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
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 lambda = context.getParameter(AmoebaVdwForce::Lambda());
if (useCutoff) {
......@@ -1090,7 +1096,6 @@ void ReferenceCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& con
throw OpenMMException("updateParametersInContext: The number of particles has changed");
// Record the values.
for (int i = 0; i < numParticles; ++i) {
int indexIV;
double sigma, epsilon, reduction;
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ 9c95f25c
......@@ -499,11 +499,11 @@ private:
int numParticles;
int useCutoff;
int usePBC;
int alchemicalMethod;
double n;
double alpha;
double cutoff;
double dispersionCoefficient;
AmoebaVdwForce::AlchemicalMethod alchemicalMethod;
int n;
double alpha;
std::vector<int> indexIVs;
std::vector< std::set<int> > allExclusions;
std::vector<double> sigmas;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp
View file @ 9c95f25c
......@@ -364,7 +364,7 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double
std::vector<Vec3> reducedPositions;
setReducedPositions(numParticles, particlePositions, indexIVs, reductions, reducedPositions);
// loop over neighbor list
// (1) calculate pair vdw ixn
// (2) accumulate forces: if particle is a site where interaction position != particle position,
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaVdwForce.cpp
View file @ 9c95f25c
......@@ -200,17 +200,26 @@ void testVdw() {
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
System system;
AmoebaVdwForce* amoebaVdwForce = new AmoebaVdwForce();;
AmoebaVdwForce* amoebaVdwForce = new AmoebaVdwForce();
int numberOfParticles = 8;
amoebaVdwForce->setSigmaCombiningRule(sigmaCombiningRule);
amoebaVdwForce->setEpsilonCombiningRule(epsilonCombiningRule);
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) {
Vec3 a(boxDimension, 0.0, 0.0);
Vec3 b(0.0, boxDimension, 0.0);
......@@ -239,16 +248,16 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co
amoebaVdwForce->addParticle(0, 1.3500000e-01, 8.3680000e-02, 9.1000000e-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);
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);
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);
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
......@@ -336,12 +345,27 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co
LangevinIntegrator integrator(0.0, 0.1, 0.01);
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);
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
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,
std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance) {
......@@ -381,6 +405,43 @@ void testVdwAmmoniaCubicMeanHhg() {
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
void testVdwAmmoniaArithmeticArithmetic() {
......@@ -2076,6 +2137,25 @@ int main(int numberOfArguments, char* argv[]) {
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) {
std::cout << "exception: " << e.what() << std::endl;
......
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