Skip to content

GitLab

Commit f2e0735f authored by Peter Eastman's avatar Peter Eastman
Browse files

Added updateParametersInContext() methods to most AMOEBA forces.

parent e6696626
Show whitespace changes
Inline Side-by-side
Showing with 1049 additions and 180 deletions
plugins/amoeba/openmmapi/src/AmoebaWcaDispersionForceImpl.cpp
View file @ f2e0735f
......@@ -141,3 +141,6 @@ std::vector<std::string> AmoebaWcaDispersionForceImpl::getKernelNames() {
return names;
}
void AmoebaWcaDispersionForceImpl::updateParametersInContext(ContextImpl& context) {
kernel.getAs<CalcAmoebaWcaDispersionForceKernel>().copyParametersToContext(context, owner);
}
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ f2e0735f
......@@ -123,6 +123,30 @@ double CudaCalcAmoebaBondForceKernel::execute(ContextImpl& context, bool include
return 0.0;
}
void CudaCalcAmoebaBondForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaBondForce& force) {
cu.setAsCurrent();
int numContexts = cu.getPlatformData().contexts.size();
int startIndex = cu.getContextIndex()*force.getNumBonds()/numContexts;
int endIndex = (cu.getContextIndex()+1)*force.getNumBonds()/numContexts;
if (numBonds != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of bonds has changed");
// Record the per-bond parameters.
vector<float2> paramVector(numBonds);
for (int i = 0; i < numBonds; i++) {
int atom1, atom2;
double length, k;
force.getBondParameters(startIndex+i, atom1, atom2, length, k);
paramVector[i] = make_float2((float) length, (float) k);
}
params->upload(paramVector);
// Mark that the current reordering may be invalid.
cu.invalidateMolecules();
}
/* -------------------------------------------------------------------------- *
* AmoebaAngleForce *
* -------------------------------------------------------------------------- */
......@@ -197,6 +221,30 @@ double CudaCalcAmoebaAngleForceKernel::execute(ContextImpl& context, bool includ
return 0.0;
}
void CudaCalcAmoebaAngleForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaAngleForce& force) {
cu.setAsCurrent();
int numContexts = cu.getPlatformData().contexts.size();
int startIndex = cu.getContextIndex()*force.getNumAngles()/numContexts;
int endIndex = (cu.getContextIndex()+1)*force.getNumAngles()/numContexts;
if (numAngles != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of angles has changed");
// Record the per-angle parameters.
vector<float2> paramVector(numAngles);
for (int i = 0; i < numAngles; i++) {
int atom1, atom2, atom3;
double angle, k;
force.getAngleParameters(startIndex+i, atom1, atom2, atom3, angle, k);
paramVector[i] = make_float2((float) angle, (float) k);
}
params->upload(paramVector);
// Mark that the current reordering may be invalid.
cu.invalidateMolecules();
}
/* -------------------------------------------------------------------------- *
* AmoebaInPlaneAngleForce *
* -------------------------------------------------------------------------- */
......@@ -271,6 +319,30 @@ double CudaCalcAmoebaInPlaneAngleForceKernel::execute(ContextImpl& context, bool
return 0.0;
}
void CudaCalcAmoebaInPlaneAngleForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaInPlaneAngleForce& force) {
cu.setAsCurrent();
int numContexts = cu.getPlatformData().contexts.size();
int startIndex = cu.getContextIndex()*force.getNumAngles()/numContexts;
int endIndex = (cu.getContextIndex()+1)*force.getNumAngles()/numContexts;
if (numAngles != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of in-plane angles has changed");
// Record the per-angle parameters.
vector<float2> paramVector(numAngles);
for (int i = 0; i < numAngles; i++) {
int atom1, atom2, atom3, atom4;
double angle, k;
force.getAngleParameters(startIndex+i, atom1, atom2, atom3, atom4, angle, k);
paramVector[i] = make_float2((float) angle, (float) k);
}
params->upload(paramVector);
// Mark that the current reordering may be invalid.
cu.invalidateMolecules();
}
/* -------------------------------------------------------------------------- *
* AmoebaPiTorsion *
* -------------------------------------------------------------------------- */
......@@ -342,6 +414,30 @@ double CudaCalcAmoebaPiTorsionForceKernel::execute(ContextImpl& context, bool in
return 0.0;
}
void CudaCalcAmoebaPiTorsionForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaPiTorsionForce& force) {
cu.setAsCurrent();
int numContexts = cu.getPlatformData().contexts.size();
int startIndex = cu.getContextIndex()*force.getNumPiTorsions()/numContexts;
int endIndex = (cu.getContextIndex()+1)*force.getNumPiTorsions()/numContexts;
if (numPiTorsions != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
// Record the per-torsion parameters.
vector<float> paramVector(numPiTorsions);
for (int i = 0; i < numPiTorsions; i++) {
int atom1, atom2, atom3, atom4, atom5, atom6;
double k;
force.getPiTorsionParameters(startIndex+i, atom1, atom2, atom3, atom4, atom5, atom6, k);
paramVector[i] = (float) k;
}
params->upload(paramVector);
// Mark that the current reordering may be invalid.
cu.invalidateMolecules();
}
/* -------------------------------------------------------------------------- *
* AmoebaStretchBend *
* -------------------------------------------------------------------------- */
......@@ -411,6 +507,30 @@ double CudaCalcAmoebaStretchBendForceKernel::execute(ContextImpl& context, bool
return 0.0;
}
void CudaCalcAmoebaStretchBendForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaStretchBendForce& force) {
cu.setAsCurrent();
int numContexts = cu.getPlatformData().contexts.size();
int startIndex = cu.getContextIndex()*force.getNumStretchBends()/numContexts;
int endIndex = (cu.getContextIndex()+1)*force.getNumStretchBends()/numContexts;
if (numStretchBends != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of bend-stretch terms has changed");
// Record the per-stretch-bend parameters.
vector<float4> paramVector(numStretchBends);
for (int i = 0; i < numStretchBends; i++) {
int atom1, atom2, atom3;
double lengthAB, lengthCB, angle, k;
force.getStretchBendParameters(startIndex+i, atom1, atom2, atom3, lengthAB, lengthCB, angle, k);
paramVector[i] = make_float4((float) lengthAB, (float) lengthCB, (float) angle, (float) k);
}
params->upload(paramVector);
// Mark that the current reordering may be invalid.
cu.invalidateMolecules();
}
/* -------------------------------------------------------------------------- *
* AmoebaOutOfPlaneBend *
* -------------------------------------------------------------------------- */
......@@ -485,6 +605,30 @@ double CudaCalcAmoebaOutOfPlaneBendForceKernel::execute(ContextImpl& context, bo
return 0.0;
}
void CudaCalcAmoebaOutOfPlaneBendForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaOutOfPlaneBendForce& force) {
cu.setAsCurrent();
int numContexts = cu.getPlatformData().contexts.size();
int startIndex = cu.getContextIndex()*force.getNumOutOfPlaneBends()/numContexts;
int endIndex = (cu.getContextIndex()+1)*force.getNumOutOfPlaneBends()/numContexts;
if (numOutOfPlaneBends != endIndex-startIndex)
throw OpenMMException("updateParametersInContext: The number of out-of-plane bends has changed");
// Record the per-bend parameters.
vector<float> paramVector(numOutOfPlaneBends);
for (int i = 0; i < numOutOfPlaneBends; i++) {
int atom1, atom2, atom3, atom4;
double k;
force.getOutOfPlaneBendParameters(startIndex+i, atom1, atom2, atom3, atom4, k);
paramVector[i] = (float) k;
}
params->upload(paramVector);
// Mark that the current reordering may be invalid.
cu.invalidateMolecules();
}
/* -------------------------------------------------------------------------- *
* AmoebaTorsionTorsion *
* -------------------------------------------------------------------------- */
......@@ -1558,6 +1702,61 @@ void CudaCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl&
computeSystemMultipoleMoments<float, float4, float4>(context, outputMultipoleMoments);
}
void CudaCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaMultipoleForce& force) {
// Make sure the new parameters are acceptable.
cu.setAsCurrent();
if (force.getNumMultipoles() != cu.getNumAtoms())
throw OpenMMException("updateParametersInContext: The number of multipoles has changed");
// Record the per-multipole parameters.
cu.getPosq().download(cu.getPinnedBuffer());
float4* posqf = (float4*) cu.getPinnedBuffer();
double4* posqd = (double4*) cu.getPinnedBuffer();
vector<float2> dampingAndTholeVec;
vector<float> polarizabilityVec;
vector<float> molecularDipolesVec;
vector<float> molecularQuadrupolesVec;
vector<int4> multipoleParticlesVec;
for (int i = 0; i < force.getNumMultipoles(); i++) {
double charge, thole, damping, polarity;
int axisType, atomX, atomY, atomZ;
vector<double> dipole, quadrupole;
force.getMultipoleParameters(i, charge, dipole, quadrupole, axisType, atomZ, atomX, atomY, thole, damping, polarity);
if (cu.getUseDoublePrecision())
posqd[i].w = charge;
else
posqf[i].w = (float) charge;
dampingAndTholeVec.push_back(make_float2((float) damping, (float) thole));
polarizabilityVec.push_back((float) polarity);
multipoleParticlesVec.push_back(make_int4(atomX, atomY, atomZ, axisType));
for (int j = 0; j < 3; j++)
molecularDipolesVec.push_back((float) dipole[j]);
molecularQuadrupolesVec.push_back((float) quadrupole[0]);
molecularQuadrupolesVec.push_back((float) quadrupole[1]);
molecularQuadrupolesVec.push_back((float) quadrupole[2]);
molecularQuadrupolesVec.push_back((float) quadrupole[4]);
molecularQuadrupolesVec.push_back((float) quadrupole[5]);
}
for (int i = force.getNumMultipoles(); i < cu.getPaddedNumAtoms(); i++) {
dampingAndTholeVec.push_back(make_float2(0, 0));
polarizabilityVec.push_back(0);
multipoleParticlesVec.push_back(make_int4(0, 0, 0, 0));
for (int j = 0; j < 3; j++)
molecularDipolesVec.push_back(0);
for (int j = 0; j < 5; j++)
molecularQuadrupolesVec.push_back(0);
}
dampingAndThole->upload(dampingAndTholeVec);
polarizability->upload(polarizabilityVec);
multipoleParticles->upload(multipoleParticlesVec);
molecularDipoles->upload(molecularDipolesVec);
molecularQuadrupoles->upload(molecularQuadrupolesVec);
cu.getPosq().upload(cu.getPinnedBuffer());
cu.invalidateMolecules();
}
/* -------------------------------------------------------------------------- *
* AmoebaGeneralizedKirkwood *
* -------------------------------------------------------------------------- */
......@@ -1770,6 +1969,25 @@ void CudaCalcAmoebaGeneralizedKirkwoodForceKernel::finishComputation(CudaArray&
cu.executeKernel(ediffKernel, ediffArgs, numForceThreadBlocks*ediffThreads, ediffThreads);
}
void CudaCalcAmoebaGeneralizedKirkwoodForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaGeneralizedKirkwoodForce& force) {
// Make sure the new parameters are acceptable.
cu.setAsCurrent();
if (force.getNumParticles() != cu.getNumAtoms())
throw OpenMMException("updateParametersInContext: The number of particles has changed");
// Record the per-particle parameters.
vector<float2> paramsVector(cu.getPaddedNumAtoms());
for (int i = 0; i < force.getNumParticles(); i++) {
double charge, radius, scalingFactor;
force.getParticleParameters(i, charge, radius, scalingFactor);
paramsVector[i] = make_float2((float) radius, (float) (scalingFactor*radius));
}
params->upload(paramsVector);
cu.invalidateMolecules();
}
/* -------------------------------------------------------------------------- *
* AmoebaVdw *
* -------------------------------------------------------------------------- */
......@@ -1913,6 +2131,36 @@ double CudaCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeF
return dispersionCoefficient/(box.x*box.y*box.z);
}
void CudaCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaVdwForce& force) {
// Make sure the new parameters are acceptable.
cu.setAsCurrent();
if (force.getNumParticles() != cu.getNumAtoms())
throw OpenMMException("updateParametersInContext: The number of particles has changed");
// Record the per-particle parameters.
vector<float2> sigmaEpsilonVec(cu.getPaddedNumAtoms(), make_float2(0, 1));
vector<int> bondReductionAtomsVec(cu.getPaddedNumAtoms(), 0);
vector<float> bondReductionFactorsVec(cu.getPaddedNumAtoms(), 0);
for (int i = 0; i < force.getNumParticles(); i++) {
int ivIndex;
double sigma, epsilon, reductionFactor;
force.getParticleParameters(i, ivIndex, sigma, epsilon, reductionFactor);
sigmaEpsilonVec[i] = make_float2((float) sigma, (float) epsilon);
bondReductionAtomsVec[i] = ivIndex;
bondReductionFactorsVec[i] = (float) reductionFactor;
}
sigmaEpsilon->upload(sigmaEpsilonVec);
bondReductionAtoms->upload(bondReductionAtomsVec);
bondReductionFactors->upload(bondReductionFactorsVec);
if (force.getUseDispersionCorrection())
dispersionCoefficient = AmoebaVdwForceImpl::calcDispersionCorrection(system, force);
else
dispersionCoefficient = 0.0;
cu.invalidateMolecules();
}
/* -------------------------------------------------------------------------- *
* AmoebaWcaDispersion *
* -------------------------------------------------------------------------- */
......@@ -1936,6 +2184,7 @@ CudaCalcAmoebaWcaDispersionForceKernel::CudaCalcAmoebaWcaDispersionForceKernel(s
}
CudaCalcAmoebaWcaDispersionForceKernel::~CudaCalcAmoebaWcaDispersionForceKernel() {
cu.setAsCurrent();
if (radiusEpsilon != NULL)
delete radiusEpsilon;
}
......@@ -1991,3 +2240,22 @@ double CudaCalcAmoebaWcaDispersionForceKernel::execute(ContextImpl& context, boo
cu.executeKernel(forceKernel, forceArgs, numForceThreadBlocks*forceThreadBlockSize, forceThreadBlockSize);
return totalMaximumDispersionEnergy;
}
void CudaCalcAmoebaWcaDispersionForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaWcaDispersionForce& force) {
// Make sure the new parameters are acceptable.
cu.setAsCurrent();
if (force.getNumParticles() != cu.getNumAtoms())
throw OpenMMException("updateParametersInContext: The number of particles has changed");
// Record the per-particle parameters.
vector<float2> radiusEpsilonVec(cu.getPaddedNumAtoms(), make_float2(0, 0));
for (int i = 0; i < cu.getNumAtoms(); i++) {
double radius, epsilon;
force.getParticleParameters(i, radius, epsilon);
radiusEpsilonVec[i] = make_float2((float) radius, (float) epsilon);
}
radiusEpsilon->upload(radiusEpsilonVec);
cu.invalidateMolecules();
}
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h
View file @ f2e0735f
......@@ -65,6 +65,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaBondForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaBondForce& force);
private:
class ForceInfo;
int numBonds;
......@@ -96,6 +103,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaAngleForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaAngleForce& force);
private:
class ForceInfo;
int numAngles;
......@@ -127,6 +141,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaInPlaneAngleForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaInPlaneAngleForce& force);
private:
class ForceInfo;
int numAngles;
......@@ -158,6 +179,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaPiTorsionForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaPiTorsionForce& force);
private:
class ForceInfo;
int numPiTorsions;
......@@ -189,6 +217,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaStretchBendForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaStretchBendForce& force);
private:
class ForceInfo;
int numStretchBends;
......@@ -220,6 +255,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaOutOfPlaneBendForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaOutOfPlaneBendForce& force);
private:
class ForceInfo;
int numOutOfPlaneBends;
......@@ -306,8 +348,13 @@ public:
* quadrupole_zx, quadrupole_zy, quadrupole_zz )
*/
void getSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaMultipoleForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaMultipoleForce& force);
private:
class ForceInfo;
class SortTrait : public CudaSort::SortTrait {
......@@ -419,6 +466,13 @@ public:
CudaArray* getInducedDipolesPolar() {
return inducedDipolePolarS;
}
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaGeneralizedKirkwoodForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaGeneralizedKirkwoodForce& force);
private:
class ForceInfo;
CudaContext& cu;
......@@ -460,6 +514,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaVdwForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaVdwForce& force);
private:
class ForceInfo;
CudaContext& cu;
......@@ -498,6 +559,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaWcaDispersionForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaWcaDispersionForce& force);
private:
class ForceInfo;
CudaContext& cu;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaAngleForce.cpp
View file @ f2e0735f
......@@ -300,6 +300,20 @@ void testOneAngle( FILE* log ) {
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log );
// Try changing the angle parameters and make sure it's still correct.
amoebaAngleForce->setAngleParameters(0, 0, 1, 2, 1.1*angle, 1.4*quadraticK);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log );
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaAngleForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log );
}
int main(int argc, char* argv[]) {
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaBondForce.cpp
View file @ f2e0735f
......@@ -204,6 +204,22 @@ void testTwoBond( FILE* log ) {
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log );
// Try changing the bond parameters and make sure it's still correct.
amoebaBondForce->setBondParameters(0, 0, 1, 1.1*bondLength, 1.4*quadraticK);
amoebaBondForce->setBondParameters(1, 1, 2, 1.2*bondLength, 0.9*quadraticK);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log );
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaBondForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log );
}
int main(int argc, char* argv[]) {
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp
View file @ f2e0735f
......@@ -30,7 +30,7 @@
* -------------------------------------------------------------------------- */
 
/**
* This tests the CUDA implementation of CudaAmoebaMultipoleForce.
* This tests the CUDA implementation of AmoebaGeneralizedKirkwoodForce.
*/
 
#include "openmm/internal/AssertionUtilities.h"
......@@ -56,14 +56,11 @@ extern "C" void registerAmoebaCudaKernelFactories();
 
// setup for 2 ammonia molecules
 
static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::PolarizationType polarizationType,
int includeCavityTerm, std::vector<Vec3>& forces, double& energy, FILE* log ){
static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce,
AmoebaMultipoleForce::PolarizationType polarizationType, int includeCavityTerm) {
 
// beginning of Multipole setup
 
System system;
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
int numberOfParticles = 8;
 
......@@ -261,7 +258,6 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Polar
 
// GK force
 
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01 );
amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00 );
amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm( includeCavityTerm );
......@@ -275,25 +271,10 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Polar
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 );
}
system.addForce(amoebaGeneralizedKirkwoodForce);
}
 
// 1-2 bonds needed
/*
AmoebaBondForce* AmoebaBondForce = new AmoebaBondForce();
// addBond: particle1, particle2, length, quadraticK
AmoebaBondForce->addBond( 0, 1, 0.0000000e+00, 0.0000000e+00 );
AmoebaBondForce->addBond( 0, 2, 0.0000000e+00, 0.0000000e+00 );
AmoebaBondForce->addBond( 0, 3, 0.0000000e+00, 0.0000000e+00 );
AmoebaBondForce->addBond( 4, 5, 0.0000000e+00, 0.0000000e+00 );
AmoebaBondForce->addBond( 4, 6, 0.0000000e+00, 0.0000000e+00 );
AmoebaBondForce->addBond( 4, 7, 0.0000000e+00, 0.0000000e+00 );
AmoebaBondForce->setAmoebaGlobalBondCubic( -2.5500000e+01 );
AmoebaBondForce->setAmoebaGlobalBondQuartic( 3.7931250e+02 );
system.addForce(AmoebaBondForce);
*/
std::vector<Vec3> positions(numberOfParticles);
static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& forces, double& energy, FILE* log) {
std::vector<Vec3> positions(context.getSystem().getNumParticles());
 
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03 );
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02 );
......@@ -304,11 +285,6 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Polar
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02 );
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02 );
 
std::string platformName;
platformName = "CUDA";
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName( platformName ) );
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
......@@ -7001,8 +6977,8 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
// compare forces and energies
 
static void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance, FILE* log ) {
const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) {
 
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
......@@ -7128,7 +7104,12 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) {
std::vector<Vec3> forces;
double energy;
 
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Direct, 0, forces, energy, log );
System system;
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
 
double expectedEnergy = -7.6636680e+01;
......@@ -7156,7 +7137,12 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) {
std::vector<Vec3> forces;
double energy;
 
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Mutual, 0, forces, energy, log );
System system;
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
 
double expectedEnergy = -7.8018875e+01;
......@@ -7175,6 +7161,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) {
}
 
// test GK mutual polarization for system comprised of two ammonia molecules
// including cavity term
 
static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE* log ) {
 
......@@ -7184,7 +7171,12 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE
std::vector<Vec3> forces;
double energy;
 
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Mutual, 1, forces, energy, log );
System system;
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 1);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
 
double expectedEnergy = -6.0434582e+01;
......@@ -7200,6 +7192,31 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE
 
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
// Try changing the particle parameters and make sure it's still correct.
for (int i = 0; i < numberOfParticles; i++) {
double charge, radius, scale;
amoebaGeneralizedKirkwoodForce->getParticleParameters(i, charge, radius, scale);
amoebaGeneralizedKirkwoodForce->setParticleParameters(i, charge, 0.9*radius, 1.1*scale);
}
LangevinIntegrator integrator2(0.0, 0.1, 0.01);
Context context2(system, integrator2, context.getPlatform());
context2.setPositions(context.getState(State::Positions).getPositions());
State state1 = context.getState(State::Forces | State::Energy);
State state2 = context2.getState(State::Forces | State::Energy);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log);
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaGeneralizedKirkwoodForce->updateParametersInContext(context);
state1 = context.getState(State::Forces | State::Energy);
compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log);
}
 
// test GK direct polarization for villin system
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaInPlaneAngleForce.cpp
View file @ f2e0735f
......@@ -377,6 +377,20 @@ void testOneAngle( FILE* log ) {
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log );
// Try changing the angle parameters and make sure it's still correct.
amoebaInPlaneAngleForce->setAngleParameters(0, 0, 1, 2, 3, 1.1*angle, 1.4*quadraticK);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log );
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaInPlaneAngleForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log );
}
int main(int argc, char* argv[]) {
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaMultipoleForce.cpp
View file @ f2e0735f
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
......@@ -30,7 +30,7 @@
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of CudaAmoebaMultipoleForce.
* This tests the CUDA implementation of AmoebaMultipoleForce.
*/
#include "openmm/internal/AssertionUtilities.h"
......@@ -56,13 +56,9 @@ extern "C" void registerAmoebaCudaKernelFactories();
// setup for 2 ammonia molecules
static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::NonbondedMethod nonbondedMethod,
static void setupMultipoleAmmonia(System& system, AmoebaMultipoleForce* amoebaMultipoleForce, AmoebaMultipoleForce::NonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::PolarizationType polarizationType,
double cutoff, int inputPmeGridDimension, std::vector<Vec3>& forces, double& energy, FILE* log ){
// beginning of Multipole setup
System system;
double cutoff, int inputPmeGridDimension) {
// box
......@@ -72,7 +68,6 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Nonbo
Vec3 c( 0.0, 0.0, boxDimension );
system.setDefaultPeriodicBoxVectors( a, b, c );
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
int numberOfParticles = 8;
amoebaMultipoleForce->setNonbondedMethod( nonbondedMethod );
......@@ -130,23 +125,23 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Nonbo
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 2, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 2, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addParticle( 1.932e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 2, 3, 3.9e-01, 2.8135002e-01, 4.96e-04 );
amoebaMultipoleForce->addParticle( 1.932e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 3, 3.9e-01, 2.8135002e-01, 4.96e-04 );
amoebaMultipoleForce->addParticle( 1.932e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 2, 3.9e-01, 2.8135002e-01, 4.96e-04 );
// second N
system.addParticle( 1.4007000e+01 );
amoebaMultipoleForce->addParticle( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 5, 6, 7, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03 );
amoebaMultipoleForce->addParticle( -5.796e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 5, 6, 7, 3.9e-01, 3.1996314e-01, 1.073e-03 );
// 3 H attached to second N
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 6, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 6, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addParticle( 1.932e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 6, 7, 3.9e-01, 2.8135002e-01, 4.96e-04 );
amoebaMultipoleForce->addParticle( 1.932e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 7, 3.9e-01, 2.8135002e-01, 4.96e-04 );
amoebaMultipoleForce->addParticle( 1.932e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 6, 3.9e-01, 2.8135002e-01, 4.96e-04 );
// covalent maps
......@@ -272,25 +267,11 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Nonbo
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
system.addForce(amoebaMultipoleForce);
}
// 1-2 bonds needed
AmoebaBondForce* amoebaBondForce = new AmoebaBondForce();
// addBond: particle1, particle2, length, quadraticK
amoebaBondForce->addBond( 0, 1, 0.0000000e+00, 0.0000000e+00 );
amoebaBondForce->addBond( 0, 2, 0.0000000e+00, 0.0000000e+00 );
amoebaBondForce->addBond( 0, 3, 0.0000000e+00, 0.0000000e+00 );
amoebaBondForce->addBond( 4, 5, 0.0000000e+00, 0.0000000e+00 );
amoebaBondForce->addBond( 4, 6, 0.0000000e+00, 0.0000000e+00 );
amoebaBondForce->addBond( 4, 7, 0.0000000e+00, 0.0000000e+00 );
amoebaBondForce->setAmoebaGlobalBondCubic( -2.5500000e+01 );
amoebaBondForce->setAmoebaGlobalBondQuartic( 3.7931250e+02 );
system.addForce(amoebaBondForce);
std::vector<Vec3> positions(numberOfParticles);
static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& forces, double& energy) {
std::vector<Vec3> positions(context.getSystem().getNumParticles());
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03 );
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02 );
......@@ -301,13 +282,6 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Nonbo
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02 );
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02 );
system.addForce(amoebaMultipoleForce);
std::string platformName;
platformName = "CUDA";
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName( platformName ) );
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
......@@ -317,9 +291,8 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Nonbo
// compare forces and energies
static void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance, FILE* log ) {
const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) {
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
......@@ -370,7 +343,6 @@ static void compareForceNormsEnergy( std::string& testName, double expectedEnerg
std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) {
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
if( log ){
......@@ -447,8 +419,13 @@ static void testMultipoleAmmoniaDirectPolarization( FILE* log ) {
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Direct,
cutoff, inputPmeGridDimension, forces, energy, log );
System system;
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
setupMultipoleAmmonia(system, amoebaMultipoleForce, AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Direct,
cutoff, inputPmeGridDimension);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
getForcesEnergyMultipoleAmmonia(context, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = -1.7428832e+01;
......@@ -478,8 +455,13 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) {
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Mutual,
cutoff, inputPmeGridDimension, forces, energy, log );
System system;
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
setupMultipoleAmmonia(system, amoebaMultipoleForce, AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Mutual,
cutoff, inputPmeGridDimension);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
getForcesEnergyMultipoleAmmonia(context, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = -1.7790449e+01;
......@@ -495,6 +477,38 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) {
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
// Try changing the particle parameters and make sure it's still correct.
for (int i = 0; i < numberOfParticles; i++) {
double charge, thole, damping, polarity;
int axisType, atomX, atomY, atomZ;
std::vector<double> dipole, quadrupole;
amoebaMultipoleForce->getMultipoleParameters(i, charge, dipole, quadrupole, axisType, atomZ, atomX, atomY, thole, damping, polarity);
dipole[0] *= 0.7;
quadrupole[2] *= 1.5;
quadrupole[6] *= 1.5;
amoebaMultipoleForce->setMultipoleParameters(i, 1.1*charge, dipole, quadrupole, axisType, atomZ, atomX, atomY, 1.3*thole, 1.4*damping, 1.5*polarity);
}
LangevinIntegrator integrator2(0.0, 0.1, 0.01);
Context context2(system, integrator2, context.getPlatform());
context2.setPositions(context.getState(State::Positions).getPositions());
State state1 = context.getState(State::Forces | State::Energy);
State state2 = context2.getState(State::Forces | State::Energy);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareForcesEnergy( testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log );
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], tolerance);
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaMultipoleForce->updateParametersInContext(context);
state1 = context.getState(State::Forces | State::Energy);
compareForcesEnergy( testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log );
}
// setup for box of 4 water molecules -- used to test PME
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaOutOfPlaneBendForce.cpp
View file @ f2e0735f
......@@ -336,6 +336,20 @@ void testOneOutOfPlaneBend( FILE* log ) {
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log );
// Try changing the bend parameters and make sure it's still correct.
amoebaOutOfPlaneBendForce->setOutOfPlaneBendParameters(0, 0, 1, 2, 3, 1.1*kOutOfPlaneBend);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log );
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaOutOfPlaneBendForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log );
}
void testOneOutOfPlaneBend2( FILE* log, int setId ) {
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaPiTorsionForce.cpp
View file @ f2e0735f
......@@ -303,6 +303,20 @@ void testOnePiTorsion( FILE* log ) {
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log );
// Try changing the torsion parameters and make sure it's still correct.
amoebaPiTorsionForce->setPiTorsionParameters(0, 0, 1, 2, 3, 4, 5, 1.2*kTorsion);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log );
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaPiTorsionForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log );
}
int main(int argc, char* argv[]) {
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaStretchBendForce.cpp
View file @ f2e0735f
......@@ -288,6 +288,20 @@ void testOneStretchBend( FILE* log ) {
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log );
// Try changing the stretch-bend parameters and make sure it's still correct.
amoebaStretchBendForce->setStretchBendParameters(0, 0, 1, 2, 1.1*abLength, 1.2*cbLength, 1.3*angleStretchBend, 1.4*kStretchBend);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log );
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaStretchBendForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log );
}
int main(int argc, char* argv[]) {
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaVdwForce.cpp
View file @ f2e0735f
......@@ -174,6 +174,35 @@ void testVdw( FILE* log ) {
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
}
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
// Try changing the particle parameters and make sure it's still correct.
for (int i = 0; i < numberOfParticles; i++) {
int indexIV;
double mass, sigma, epsilon, reduction;
amoebaVdwForce->getParticleParameters(i, indexIV, sigma, epsilon, reduction);
amoebaVdwForce->setParticleParameters(i, indexIV, 0.9*sigma, 2.0*epsilon, 0.95*reduction);
}
LangevinIntegrator integrator2(0.0, 0.1, 0.01);
Context context2(system, integrator2, Platform::getPlatformByName(platformName));
context2.setPositions(positions);
State state1 = context.getState(State::Forces | State::Energy);
State state2 = context2.getState(State::Forces | State::Energy);
bool exceptionThrown = false;
try {
// This should throw an exception.
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], tolerance);
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaVdwForce->updateParametersInContext(context);
state1 = context.getState(State::Forces | State::Energy);
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], tolerance);
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), tolerance);
}
void setupAndGetForcesEnergyVdwAmmonia( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff,
......
plugins/amoeba/platforms/cuda/tests/TestCudaWcaDispersionForce.cpp
View file @ f2e0735f
......@@ -54,13 +54,39 @@ const double TOL = 1e-4;
extern "C" void registerAmoebaCudaKernelFactories();
void setupAndGetForcesEnergyWcaDispersionAmmonia( std::vector<Vec3>& forces, double& energy, FILE* log ){
void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy,
const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e\n", testName.c_str(), expectedEnergy, energy );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC_MOD( expectedForces[ii], forces[ii], tolerance, testName );
}
ASSERT_EQUAL_TOL_MOD( expectedEnergy, energy, tolerance, testName );
}
// test Wca dispersion
void testWcaDispersionAmmonia( FILE* log ) {
std::string testName = "testWcaDispersionAmmonia";
// beginning of WcaDispersion setup
int numberOfParticles = 8;
// Create the system.
System system;
AmoebaWcaDispersionForce* amoebaWcaDispersionForce = new AmoebaWcaDispersionForce();;
int numberOfParticles = 8;
amoebaWcaDispersionForce->setEpso( 4.6024000e-01 );
amoebaWcaDispersionForce->setEpsh( 5.6484000e-02 );
......@@ -107,46 +133,12 @@ void setupAndGetForcesEnergyWcaDispersionAmmonia( std::vector<Vec3>& forces, dou
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
energy = state.getPotentialEnergy();
}
void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance, FILE* log ) {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e\n", testName.c_str(), expectedEnergy, energy );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC_MOD( expectedForces[ii], forces[ii], tolerance, testName );
}
ASSERT_EQUAL_TOL_MOD( expectedEnergy, energy, tolerance, testName );
}
// test Wca dispersion
void testWcaDispersionAmmonia( FILE* log ) {
std::string testName = "testWcaDispersionAmmonia";
int numberOfParticles = 8;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyWcaDispersionAmmonia( forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
std::vector<Vec3> forces = state.getForces();
double energy = state.getPotentialEnergy();
// TINKER-computed values
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = -2.6981209e+01;
expectedForces[0] = Vec3( 4.7839388e+00, -7.3510133e-04, -5.0382764e-01 );
......@@ -160,6 +152,31 @@ void testWcaDispersionAmmonia( FILE* log ) {
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
// Try changing the particle parameters and make sure it's still correct.
for (int i = 0; i < numberOfParticles; i++) {
double radius, epsilon;
amoebaWcaDispersionForce->getParticleParameters(i, radius, epsilon);
amoebaWcaDispersionForce->setParticleParameters(i, 0.9*radius, 2.0*epsilon);
}
LangevinIntegrator integrator2(0.0, 0.1, 0.01);
Context context2(system, integrator2, Platform::getPlatformByName(platformName));
context2.setPositions(positions);
State state1 = context.getState(State::Forces | State::Energy);
State state2 = context2.getState(State::Forces | State::Energy);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareForcesEnergy(testName, state1.getPotentialEnergy(), state2.getPotentialEnergy(), state1.getForces(), state2.getForces(), tolerance, log);
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaWcaDispersionForce->updateParametersInContext(context);
state1 = context.getState(State::Forces | State::Energy);
compareForcesEnergy(testName, state1.getPotentialEnergy(), state2.getPotentialEnergy(), state1.getForces(), state2.getForces(), tolerance, log);
}
int main(int argc, char* argv[]) {
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ f2e0735f
......@@ -111,6 +111,23 @@ double ReferenceCalcAmoebaBondForceKernel::execute(ContextImpl& context, bool in
return static_cast<double>(energy);
}
void ReferenceCalcAmoebaBondForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaBondForce& force) {
if (numBonds != force.getNumBonds())
throw OpenMMException("updateParametersInContext: The number of bonds has changed");
// Record the values.
for (int i = 0; i < numBonds; ++i) {
int particle1Index, particle2Index;
double lengthValue, kValue;
force.getBondParameters(i, particle1Index, particle2Index, lengthValue, kValue);
if (particle1Index != particle1[i] || particle2Index != particle2[i])
throw OpenMMException("updateParametersInContext: The set of particles in a bond has changed");
length[i] = (RealOpenMM) lengthValue;
kQuadratic[i] = (RealOpenMM) kValue;
}
}
// ***************************************************************************
ReferenceCalcAmoebaAngleForceKernel::ReferenceCalcAmoebaAngleForceKernel(std::string name, const Platform& platform, System& system) :
......@@ -149,6 +166,23 @@ double ReferenceCalcAmoebaAngleForceKernel::execute(ContextImpl& context, bool i
return static_cast<double>(energy);
}
void ReferenceCalcAmoebaAngleForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaAngleForce& force) {
if (numAngles != force.getNumAngles())
throw OpenMMException("updateParametersInContext: The number of angles has changed");
// Record the values.
for (int i = 0; i < numAngles; ++i) {
int particle1Index, particle2Index, particle3Index;
double angleValue, k;
force.getAngleParameters(i, particle1Index, particle2Index, particle3Index, angleValue, k);
if (particle1Index != particle1[i] || particle2Index != particle2[i] || particle3Index != particle3[i])
throw OpenMMException("updateParametersInContext: The set of particles in an angle has changed");
angle[i] = (RealOpenMM) angleValue;
kQuadratic[i] = (RealOpenMM) k;
}
}
ReferenceCalcAmoebaInPlaneAngleForceKernel::ReferenceCalcAmoebaInPlaneAngleForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaInPlaneAngleForceKernel(name, platform), system(system) {
}
......@@ -187,6 +221,23 @@ double ReferenceCalcAmoebaInPlaneAngleForceKernel::execute(ContextImpl& context,
return static_cast<double>(energy);
}
void ReferenceCalcAmoebaInPlaneAngleForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaInPlaneAngleForce& force) {
if (numAngles != force.getNumAngles())
throw OpenMMException("updateParametersInContext: The number of angles has changed");
// Record the values.
for (int i = 0; i < numAngles; ++i) {
int particle1Index, particle2Index, particle3Index, particle4Index;
double angleValue, k;
force.getAngleParameters(i, particle1Index, particle2Index, particle3Index, particle4Index, angleValue, k);
if (particle1Index != particle1[i] || particle2Index != particle2[i] || particle3Index != particle3[i] || particle4Index != particle4[i])
throw OpenMMException("updateParametersInContext: The set of particles in an angle has changed");
angle[i] = (RealOpenMM) angleValue;
kQuadratic[i] = (RealOpenMM) k;
}
}
ReferenceCalcAmoebaPiTorsionForceKernel::ReferenceCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaPiTorsionForceKernel(name, platform), system(system) {
}
......@@ -222,6 +273,23 @@ double ReferenceCalcAmoebaPiTorsionForceKernel::execute(ContextImpl& context, bo
return static_cast<double>(energy);
}
void ReferenceCalcAmoebaPiTorsionForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaPiTorsionForce& force) {
if (numPiTorsions != force.getNumPiTorsions())
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
// Record the values.
for (int i = 0; i < numPiTorsions; ++i) {
int particle1Index, particle2Index, particle3Index, particle4Index, particle5Index, particle6Index;
double kTorsionParameter;
force.getPiTorsionParameters(i, particle1Index, particle2Index, particle3Index, particle4Index, particle5Index, particle6Index, kTorsionParameter);
if (particle1Index != particle1[i] || particle2Index != particle2[i] || particle3Index != particle3[i] ||
particle4Index != particle4[i] || particle5Index != particle5[i] || particle6Index != particle6[i])
throw OpenMMException("updateParametersInContext: The set of particles in a torsion has changed");
kTorsion[i] = (RealOpenMM) kTorsionParameter;
}
}
ReferenceCalcAmoebaStretchBendForceKernel::ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaStretchBendForceKernel(name, platform), system(system) {
}
......@@ -255,6 +323,25 @@ double ReferenceCalcAmoebaStretchBendForceKernel::execute(ContextImpl& context,
return static_cast<double>(energy);
}
void ReferenceCalcAmoebaStretchBendForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaStretchBendForce& force) {
if (numStretchBends != force.getNumStretchBends())
throw OpenMMException("updateParametersInContext: The number of stretch-bends has changed");
// Record the values.
for (int i = 0; i < numStretchBends; ++i) {
int particle1Index, particle2Index, particle3Index;
double lengthAB, lengthCB, angle, k;
force.getStretchBendParameters(i, particle1Index, particle2Index, particle3Index, lengthAB, lengthCB, angle, k);
if (particle1Index != particle1[i] || particle2Index != particle2[i] || particle3Index != particle3[i])
throw OpenMMException("updateParametersInContext: The set of particles in a stretch-bend has changed");
lengthABParameters[i] = (RealOpenMM) lengthAB;
lengthCBParameters[i] = (RealOpenMM) lengthCB;
angleParameters[i] = (RealOpenMM) angle;
kParameters[i] = (RealOpenMM) k;
}
}
ReferenceCalcAmoebaOutOfPlaneBendForceKernel::ReferenceCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaOutOfPlaneBendForceKernel(name, platform), system(system) {
}
......@@ -298,6 +385,22 @@ double ReferenceCalcAmoebaOutOfPlaneBendForceKernel::execute(ContextImpl& contex
return static_cast<double>(energy);
}
void ReferenceCalcAmoebaOutOfPlaneBendForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaOutOfPlaneBendForce& force) {
if (numOutOfPlaneBends != force.getNumOutOfPlaneBends())
throw OpenMMException("updateParametersInContext: The number of out-of-plane bends has changed");
// Record the values.
for (int i = 0; i < numOutOfPlaneBends; ++i) {
int particle1Index, particle2Index, particle3Index, particle4Index;
double k;
force.getOutOfPlaneBendParameters(i, particle1Index, particle2Index, particle3Index, particle4Index, k);
if (particle1Index != particle1[i] || particle2Index != particle2[i] || particle3Index != particle3[i] || particle4Index != particle4[i])
throw OpenMMException("updateParametersInContext: The set of particles in an out-of-plane bend has changed");
kParameters[i] = (RealOpenMM) k;
}
}
ReferenceCalcAmoebaTorsionTorsionForceKernel::ReferenceCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaTorsionTorsionForceKernel(name, platform), system(system) {
}
......@@ -627,6 +730,43 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextI
return;
}
void ReferenceCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaMultipoleForce& force) {
if (numMultipoles != force.getNumMultipoles())
throw OpenMMException("updateParametersInContext: The number of multipoles has changed");
// Record the values.
int dipoleIndex = 0;
int quadrupoleIndex = 0;
for (int i = 0; i < numMultipoles; ++i) {
int axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY;
double charge, tholeD, dampingFactorD, polarityD;
std::vector<double> dipolesD;
std::vector<double> quadrupolesD;
force.getMultipoleParameters(i, charge, dipolesD, quadrupolesD, axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY, tholeD, dampingFactorD, polarityD);
axisTypes[i] = axisType;
multipoleAtomZs[i] = multipoleAtomZ;
multipoleAtomXs[i] = multipoleAtomX;
multipoleAtomYs[i] = multipoleAtomY;
charges[i] = (RealOpenMM) charge;
tholes[i] = (RealOpenMM) tholeD;
dampingFactors[i] = (RealOpenMM) dampingFactorD;
polarity[i] = (RealOpenMM) polarityD;
dipoles[dipoleIndex++] = (RealOpenMM) dipolesD[0];
dipoles[dipoleIndex++] = (RealOpenMM) dipolesD[1];
dipoles[dipoleIndex++] = (RealOpenMM) dipolesD[2];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[0];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[1];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[2];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[3];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[4];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[5];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[6];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[7];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[8];
}
}
/* -------------------------------------------------------------------------- *
* AmoebaGeneralizedKirkwood *
* -------------------------------------------------------------------------- */
......@@ -737,6 +877,21 @@ double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::execute(ContextImpl& c
return 0.0;
}
void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaGeneralizedKirkwoodForce& force) {
if (numParticles != force.getNumParticles())
throw OpenMMException("updateParametersInContext: The number of particles has changed");
// Record the values.
for (int i = 0; i < numParticles; ++i) {
double particleCharge, particleRadius, scalingFactor;
force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor);
atomicRadii[i] = particleRadius;
scaleFactors[i] = scalingFactor;
charges[i] = particleCharge;
}
}
ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaVdwForceKernel(name, platform), system(system) {
useCutoff = 0;
......@@ -819,6 +974,23 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
return static_cast<double>(energy);
}
void ReferenceCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaVdwForce& force) {
if (numParticles != force.getNumParticles())
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;
force.getParticleParameters(i, indexIV, sigma, epsilon, reduction);
indexIVs[i] = indexIV;
sigmas[i] = (RealOpenMM) sigma;
epsilons[i] = (RealOpenMM) epsilon;
reductions[i]= (RealOpenMM) reduction;
}
}
/* -------------------------------------------------------------------------- *
* AmoebaWcaDispersion *
* -------------------------------------------------------------------------- */
......@@ -865,3 +1037,17 @@ double ReferenceCalcAmoebaWcaDispersionForceKernel::execute(ContextImpl& context
RealOpenMM energy = amoebaReferenceWcaDispersionForce.calculateForceAndEnergy( numParticles, posData, radii, epsilons, totalMaximumDispersionEnergy, forceData);
return static_cast<double>(energy);
}
void ReferenceCalcAmoebaWcaDispersionForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaWcaDispersionForce& force) {
if (numParticles != force.getNumParticles())
throw OpenMMException("updateParametersInContext: The number of particles has changed");
// Record the values.
for (int i = 0; i < numParticles; ++i) {
double radius, epsilon;
force.getParticleParameters(i, radius, epsilon);
radii[i] = (RealOpenMM) radius;
epsilons[i] = (RealOpenMM) epsilon;
}
}
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ f2e0735f
......@@ -61,6 +61,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaBondForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaBondForce& force);
private:
int numBonds;
std::vector<int> particle1;
......@@ -95,6 +102,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaAngleForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaAngleForce& force);
private:
int numAngles;
std::vector<int> particle1;
......@@ -132,6 +146,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaInPlaneAngleForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaInPlaneAngleForce& force);
private:
int numAngles;
std::vector<int> particle1;
......@@ -170,6 +191,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaPiTorsionForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaPiTorsionForce& force);
private:
int numPiTorsions;
std::vector<int> particle1;
......@@ -205,6 +233,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaStretchBendForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaStretchBendForce& force);
private:
int numStretchBends;
std::vector<int> particle1;
......@@ -240,6 +275,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaOutOfPlaneBendForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaOutOfPlaneBendForce& force);
private:
int numOutOfPlaneBends;
std::vector<int> particle1;
......@@ -346,6 +388,13 @@ public:
quadrupole_zx, quadrupole_zy, quadrupole_zz )
*/
void getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMoments);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaMultipoleForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaMultipoleForce& force);
private:
......@@ -398,6 +447,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaVdwForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaVdwForce& force);
private:
int numParticles;
int useCutoff;
......@@ -438,6 +494,13 @@ public:
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaWcaDispersionForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaWcaDispersionForce& force);
private:
int numParticles;
......@@ -565,6 +628,14 @@ public:
*/
void getCharges( std::vector<RealOpenMM>& charges ) const;
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the AmoebaGeneralizedKirkwoodForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const AmoebaGeneralizedKirkwoodForce& force);
private:
int numParticles;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaAngleForce.cpp
View file @ f2e0735f
......@@ -295,6 +295,20 @@ void testOneAngle( FILE* log ) {
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log );
// Try changing the angle parameters and make sure it's still correct.
amoebaAngleForce->setAngleParameters(0, 0, 1, 2, 1.1*angle, 1.4*quadraticK);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log );
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaAngleForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log );
}
int main( int numberOfArguments, char* argv[] ) {
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaBondForce.cpp
View file @ f2e0735f
......@@ -40,6 +40,7 @@
#include "openmm/LangevinIntegrator.h"
#include <iostream>
#include <vector>
#include <math.h>
using namespace OpenMM;
......@@ -199,6 +200,22 @@ void testTwoBond( FILE* log ) {
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log );
// Try changing the bond parameters and make sure it's still correct.
amoebaBondForce->setBondParameters(0, 0, 1, 1.1*bondLength, 1.4*quadraticK);
amoebaBondForce->setBondParameters(1, 1, 2, 1.2*bondLength, 0.9*quadraticK);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log );
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaBondForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log );
}
int main( int numberOfArguments, char* argv[] ) {
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaGeneralizedKirkwoodForce.cpp
View file @ f2e0735f
......@@ -54,13 +54,11 @@ const double TOL = 1e-4;
 
// setup for 2 ammonia molecules
 
static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::PolarizationType polarizationType,
int includeCavityTerm, std::vector<Vec3>& forces, double& energy, FILE* log ){
static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce,
AmoebaMultipoleForce::PolarizationType polarizationType, int includeCavityTerm) {
 
// beginning of Multipole setup
 
System system;
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
int numberOfParticles = 8;
 
......@@ -259,7 +257,6 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Polar
 
// GK force
 
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01 );
amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00 );
amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm( includeCavityTerm );
......@@ -273,8 +270,10 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Polar
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 );
}
system.addForce(amoebaGeneralizedKirkwoodForce);
}
 
std::vector<Vec3> positions(numberOfParticles);
static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& forces, double& energy, FILE* log) {
std::vector<Vec3> positions(context.getSystem().getNumParticles());
 
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03 );
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02 );
......@@ -285,11 +284,6 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Polar
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02 );
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02 );
 
std::string platformName;
platformName = "Reference";
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName( platformName ) );
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
......@@ -6982,8 +6976,8 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
// compare forces and energies
 
static void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance, FILE* log ) {
const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) {
 
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
......@@ -7109,7 +7103,12 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) {
std::vector<Vec3> forces;
double energy;
 
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Direct, 0, forces, energy, log );
System system;
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
 
double expectedEnergy = -7.6636680e+01;
......@@ -7137,7 +7136,12 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) {
std::vector<Vec3> forces;
double energy;
 
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Mutual, 0, forces, energy, log );
System system;
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
 
double expectedEnergy = -7.8018875e+01;
......@@ -7166,7 +7170,12 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE
std::vector<Vec3> forces;
double energy;
 
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Mutual, 1, forces, energy, log );
System system;
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 1);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
 
double expectedEnergy = -6.0434582e+01;
......@@ -7182,6 +7191,31 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE
 
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
// Try changing the particle parameters and make sure it's still correct.
for (int i = 0; i < numberOfParticles; i++) {
double charge, radius, scale;
amoebaGeneralizedKirkwoodForce->getParticleParameters(i, charge, radius, scale);
amoebaGeneralizedKirkwoodForce->setParticleParameters(i, charge, 0.9*radius, 1.1*scale);
}
LangevinIntegrator integrator2(0.0, 0.1, 0.01);
Context context2(system, integrator2, context.getPlatform());
context2.setPositions(context.getState(State::Positions).getPositions());
State state1 = context.getState(State::Forces | State::Energy);
State state2 = context2.getState(State::Forces | State::Energy);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log);
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaGeneralizedKirkwoodForce->updateParametersInContext(context);
state1 = context.getState(State::Forces | State::Energy);
compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log);
}
 
// test GK direct polarization for villin system
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaInPlaneAngleForce.cpp
View file @ f2e0735f
......@@ -371,6 +371,20 @@ void testOneAngle( FILE* log ) {
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log );
// Try changing the angle parameters and make sure it's still correct.
amoebaInPlaneAngleForce->setAngleParameters(0, 0, 1, 2, 3, 1.1*angle, 1.4*quadraticK);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log );
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaInPlaneAngleForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log );
}
int main( int numberOfArguments, char* argv[] ) {
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaMultipoleForce.cpp
View file @ f2e0735f
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -54,13 +54,9 @@ const double TOL = 1e-4;
// setup for 2 ammonia molecules
static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::NonbondedMethod nonbondedMethod,
static void setupMultipoleAmmonia(System& system, AmoebaMultipoleForce* amoebaMultipoleForce, AmoebaMultipoleForce::NonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::PolarizationType polarizationType,
double cutoff, int inputPmeGridDimension, std::vector<Vec3>& forces, double& energy, FILE* log ){
// beginning of Multipole setup
System system;
double cutoff, int inputPmeGridDimension) {
// box
......@@ -70,7 +66,6 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Nonbo
Vec3 c( 0.0, 0.0, boxDimension );
system.setDefaultPeriodicBoxVectors( a, b, c );
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
int numberOfParticles = 8;
amoebaMultipoleForce->setNonbondedMethod( nonbondedMethod );
......@@ -270,8 +265,11 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Nonbo
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
system.addForce(amoebaMultipoleForce);
}
std::vector<Vec3> positions(numberOfParticles);
static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& forces, double& energy) {
std::vector<Vec3> positions(context.getSystem().getNumParticles());
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03 );
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02 );
......@@ -282,13 +280,6 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Nonbo
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02 );
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02 );
system.addForce(amoebaMultipoleForce);
std::string platformName;
platformName = "Reference";
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName( platformName ) );
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
......@@ -298,8 +289,8 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Nonbo
// compare forces and energies
static void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance, FILE* log ) {
const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) {
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
......@@ -426,8 +417,13 @@ static void testMultipoleAmmoniaDirectPolarization( FILE* log ) {
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Direct,
cutoff, inputPmeGridDimension, forces, energy, log );
System system;
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
setupMultipoleAmmonia(system, amoebaMultipoleForce, AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Direct,
cutoff, inputPmeGridDimension);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
getForcesEnergyMultipoleAmmonia(context, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = -1.7428832e+01;
......@@ -457,8 +453,13 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) {
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Mutual,
cutoff, inputPmeGridDimension, forces, energy, log );
System system;
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
setupMultipoleAmmonia(system, amoebaMultipoleForce, AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Mutual,
cutoff, inputPmeGridDimension);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
getForcesEnergyMultipoleAmmonia(context, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = -1.7790449e+01;
......@@ -474,6 +475,36 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) {
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
// Try changing the particle parameters and make sure it's still correct.
for (int i = 0; i < numberOfParticles; i++) {
double charge, thole, damping, polarity;
int axisType, atomX, atomY, atomZ;
std::vector<double> dipole, quadrupole;
amoebaMultipoleForce->getMultipoleParameters(i, charge, dipole, quadrupole, axisType, atomZ, atomX, atomY, thole, damping, polarity);
dipole[0] *= 0.7;
quadrupole[2] *= 1.5;
quadrupole[6] *= 1.5;
amoebaMultipoleForce->setMultipoleParameters(i, 1.1*charge, dipole, quadrupole, axisType, atomZ, atomX, atomY, 1.3*thole, 1.4*damping, 1.5*polarity);
}
LangevinIntegrator integrator2(0.0, 0.1, 0.01);
Context context2(system, integrator2, context.getPlatform());
context2.setPositions(context.getState(State::Positions).getPositions());
State state1 = context.getState(State::Forces | State::Energy);
State state2 = context2.getState(State::Forces | State::Energy);
bool exceptionThrown = false;
try {
// This should throw an exception.
compareForcesEnergy( testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log );
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaMultipoleForce->updateParametersInContext(context);
state1 = context.getState(State::Forces | State::Energy);
compareForcesEnergy( testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log );
}
// setup for box of 4 water molecules -- used to test PME
......
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