Skip to content

GitLab

Commit 4bc723ab authored by Peter Eastman's avatar Peter Eastman
Browse files

Implemented updateParametersInContext() for seven more Force classes

parent b5e2a951
Show whitespace changes
Inline Side-by-side
Showing with 438 additions and 68 deletions
platforms/opencl/tests/TestOpenCLPeriodicTorsionForce.cpp
View file @ 4bc723ab
......@@ -67,12 +67,29 @@ void testPeriodicTorsions() {
positions[3] = Vec3(1, 0, 2);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double torque = -2*1.1*std::sin(2*PI_M/3);
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), TOL);
}
// Try changing the torsion parameters and make sure it's still correct.
forceField->setTorsionParameters(0, 0, 1, 2, 3, 3, PI_M/3.2, 1.3);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double dtheta = (3*PI_M/2)-(PI_M/3.2);
double torque = -3*1.3*std::sin(dtheta);
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(1.3*(1+std::cos(dtheta)), state.getPotentialEnergy(), TOL);
}
}
void testParallelComputation() {
......
platforms/opencl/tests/TestOpenCLRBTorsionForce.cpp
View file @ 4bc723ab
......@@ -67,6 +67,7 @@ void testRBTorsions() {
positions[3] = Vec3(1, 1, 1);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double psi = 0.25*PI_M - PI_M;
double torque = 0.0;
......@@ -83,6 +84,31 @@ void testRBTorsions() {
energy += c*std::pow(std::cos(psi), i);
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
// Try changing the torsion parameters and make sure it's still correct.
forceField->setTorsionParameters(0, 0, 1, 2, 3, 0.11, 0.22, 0.33, 0.44, 0.55, 0.66);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double psi = 0.25*PI_M - PI_M;
double torque = 0.0;
for (int i = 1; i < 6; ++i) {
double c = 0.11*(i+1);
torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
}
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, -0.5*torque), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
double energy = 0.0;
for (int i = 0; i < 6; ++i) {
double c = 0.11*(i+1);
energy += c*std::pow(std::cos(psi), i);
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
}
void testParallelComputation() {
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 4bc723ab
......@@ -326,7 +326,7 @@ void ReferenceCalcHarmonicBondForceKernel::initialize(const System& system, cons
numBonds = force.getNumBonds();
bondIndexArray = allocateIntArray(numBonds, 2);
bondParamArray = allocateRealArray(numBonds, 2);
for (int i = 0; i < force.getNumBonds(); ++i) {
for (int i = 0; i < numBonds; ++i) {
int particle1, particle2;
double length, k;
force.getBondParameters(i, particle1, particle2, length, k);
......@@ -347,6 +347,25 @@ double ReferenceCalcHarmonicBondForceKernel::execute(ContextImpl& context, bool
return energy;
}
void ReferenceCalcHarmonicBondForceKernel::copyParametersToContext(ContextImpl& context, const HarmonicBondForce& 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 particle1, particle2;
double length, k;
force.getBondParameters(i, particle1, particle2, length, k);
if (particle1 != bondIndexArray[i][0] || particle2 != bondIndexArray[i][1])
throw OpenMMException("updateParametersInContext: The set of particles in a bond has changed");
bondIndexArray[i][0] = particle1;
bondIndexArray[i][1] = particle2;
bondParamArray[i][0] = (RealOpenMM) length;
bondParamArray[i][1] = (RealOpenMM) k;
}
}
ReferenceCalcCustomBondForceKernel::~ReferenceCalcCustomBondForceKernel() {
disposeIntArray(bondIndexArray, numBonds);
disposeRealArray(bondParamArray, numBonds);
......@@ -361,7 +380,7 @@ void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const
bondIndexArray = allocateIntArray(numBonds, 2);
bondParamArray = allocateRealArray(numBonds, numParameters);
vector<double> params;
for (int i = 0; i < force.getNumBonds(); ++i) {
for (int i = 0; i < numBonds; ++i) {
int particle1, particle2;
force.getBondParameters(i, particle1, particle2, params);
bondIndexArray[i][0] = particle1;
......@@ -394,6 +413,24 @@ double ReferenceCalcCustomBondForceKernel::execute(ContextImpl& context, bool in
return energy;
}
void ReferenceCalcCustomBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomBondForce& force) {
if (numBonds != force.getNumBonds())
throw OpenMMException("updateParametersInContext: The number of bonds has changed");
// Record the values.
int numParameters = force.getNumPerBondParameters();
vector<double> params;
for (int i = 0; i < numBonds; ++i) {
int particle1, particle2;
force.getBondParameters(i, particle1, particle2, params);
if (particle1 != bondIndexArray[i][0] || particle2 != bondIndexArray[i][1])
throw OpenMMException("updateParametersInContext: The set of particles in a bond has changed");
for (int j = 0; j < numParameters; j++)
bondParamArray[i][j] = (RealOpenMM) params[j];
}
}
ReferenceCalcHarmonicAngleForceKernel::~ReferenceCalcHarmonicAngleForceKernel() {
disposeIntArray(angleIndexArray, numAngles);
disposeRealArray(angleParamArray, numAngles);
......@@ -403,7 +440,7 @@ void ReferenceCalcHarmonicAngleForceKernel::initialize(const System& system, con
numAngles = force.getNumAngles();
angleIndexArray = allocateIntArray(numAngles, 3);
angleParamArray = allocateRealArray(numAngles, 2);
for (int i = 0; i < force.getNumAngles(); ++i) {
for (int i = 0; i < numAngles; ++i) {
int particle1, particle2, particle3;
double angle, k;
force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
......@@ -425,6 +462,23 @@ double ReferenceCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool
return energy;
}
void ReferenceCalcHarmonicAngleForceKernel::copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& 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 particle1, particle2, particle3;
double angle, k;
force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
if (particle1 != angleIndexArray[i][0] || particle2 != angleIndexArray[i][1] || particle3 != angleIndexArray[i][2])
throw OpenMMException("updateParametersInContext: The set of particles in an angle has changed");
angleParamArray[i][0] = (RealOpenMM) angle;
angleParamArray[i][1] = (RealOpenMM) k;
}
}
ReferenceCalcCustomAngleForceKernel::~ReferenceCalcCustomAngleForceKernel() {
disposeIntArray(angleIndexArray, numAngles);
disposeRealArray(angleParamArray, numAngles);
......@@ -439,7 +493,7 @@ void ReferenceCalcCustomAngleForceKernel::initialize(const System& system, const
angleIndexArray = allocateIntArray(numAngles, 3);
angleParamArray = allocateRealArray(numAngles, numParameters);
vector<double> params;
for (int i = 0; i < force.getNumAngles(); ++i) {
for (int i = 0; i < numAngles; ++i) {
int particle1, particle2, particle3;
force.getAngleParameters(i, particle1, particle2, particle3, params);
angleIndexArray[i][0] = particle1;
......@@ -473,6 +527,24 @@ double ReferenceCalcCustomAngleForceKernel::execute(ContextImpl& context, bool i
return energy;
}
void ReferenceCalcCustomAngleForceKernel::copyParametersToContext(ContextImpl& context, const CustomAngleForce& force) {
if (numAngles != force.getNumAngles())
throw OpenMMException("updateParametersInContext: The number of angles has changed");
// Record the values.
int numParameters = force.getNumPerAngleParameters();
vector<double> params;
for (int i = 0; i < numAngles; ++i) {
int particle1, particle2, particle3;
force.getAngleParameters(i, particle1, particle2, particle3, params);
if (particle1 != angleIndexArray[i][0] || particle2 != angleIndexArray[i][1] || particle3 != angleIndexArray[i][2])
throw OpenMMException("updateParametersInContext: The set of particles in an angle has changed");
for (int j = 0; j < numParameters; j++)
angleParamArray[i][j] = (RealOpenMM) params[j];
}
}
ReferenceCalcPeriodicTorsionForceKernel::~ReferenceCalcPeriodicTorsionForceKernel() {
disposeIntArray(torsionIndexArray, numTorsions);
disposeRealArray(torsionParamArray, numTorsions);
......@@ -482,7 +554,7 @@ void ReferenceCalcPeriodicTorsionForceKernel::initialize(const System& system, c
numTorsions = force.getNumTorsions();
torsionIndexArray = allocateIntArray(numTorsions, 4);
torsionParamArray = allocateRealArray(numTorsions, 3);
for (int i = 0; i < force.getNumTorsions(); ++i) {
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4, periodicity;
double phase, k;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, periodicity, phase, k);
......@@ -506,6 +578,24 @@ double ReferenceCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bo
return energy;
}
void ReferenceCalcPeriodicTorsionForceKernel::copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force) {
if (numTorsions != force.getNumTorsions())
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
// Record the values.
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4, periodicity;
double phase, k;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, periodicity, phase, k);
if (particle1 != torsionIndexArray[i][0] || particle2 != torsionIndexArray[i][1] || particle3 != torsionIndexArray[i][2] || particle4 != torsionIndexArray[i][3])
throw OpenMMException("updateParametersInContext: The set of particles in a torsion has changed");
torsionParamArray[i][0] = (RealOpenMM) k;
torsionParamArray[i][1] = (RealOpenMM) phase;
torsionParamArray[i][2] = (RealOpenMM) periodicity;
}
}
ReferenceCalcRBTorsionForceKernel::~ReferenceCalcRBTorsionForceKernel() {
disposeIntArray(torsionIndexArray, numTorsions);
disposeRealArray(torsionParamArray, numTorsions);
......@@ -515,7 +605,7 @@ void ReferenceCalcRBTorsionForceKernel::initialize(const System& system, const R
numTorsions = force.getNumTorsions();
torsionIndexArray = allocateIntArray(numTorsions, 4);
torsionParamArray = allocateRealArray(numTorsions, 6);
for (int i = 0; i < force.getNumTorsions(); ++i) {
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4;
double c0, c1, c2, c3, c4, c5;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
......@@ -542,6 +632,27 @@ double ReferenceCalcRBTorsionForceKernel::execute(ContextImpl& context, bool inc
return energy;
}
void ReferenceCalcRBTorsionForceKernel::copyParametersToContext(ContextImpl& context, const RBTorsionForce& force) {
if (numTorsions != force.getNumTorsions())
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
// Record the values.
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4;
double c0, c1, c2, c3, c4, c5;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
if (particle1 != torsionIndexArray[i][0] || particle2 != torsionIndexArray[i][1] || particle3 != torsionIndexArray[i][2] || particle4 != torsionIndexArray[i][3])
throw OpenMMException("updateParametersInContext: The set of particles in a torsion has changed");
torsionParamArray[i][0] = (RealOpenMM) c0;
torsionParamArray[i][1] = (RealOpenMM) c1;
torsionParamArray[i][2] = (RealOpenMM) c2;
torsionParamArray[i][3] = (RealOpenMM) c3;
torsionParamArray[i][4] = (RealOpenMM) c4;
torsionParamArray[i][5] = (RealOpenMM) c5;
}
}
void ReferenceCalcCMAPTorsionForceKernel::initialize(const System& system, const CMAPTorsionForce& force) {
int numMaps = force.getNumMaps();
int numTorsions = force.getNumTorsions();
......@@ -591,7 +702,7 @@ void ReferenceCalcCustomTorsionForceKernel::initialize(const System& system, con
torsionIndexArray = allocateIntArray(numTorsions, 4);
torsionParamArray = allocateRealArray(numTorsions, numParameters);
vector<double> params;
for (int i = 0; i < force.getNumTorsions(); ++i) {
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, params);
torsionIndexArray[i][0] = particle1;
......@@ -626,6 +737,24 @@ double ReferenceCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool
return energy;
}
void ReferenceCalcCustomTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CustomTorsionForce& force) {
if (numTorsions != force.getNumTorsions())
throw OpenMMException("updateParametersInContext: The number of torsions has changed");
// Record the values.
int numParameters = force.getNumPerTorsionParameters();
vector<double> params;
for (int i = 0; i < numTorsions; ++i) {
int particle1, particle2, particle3, particle4;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, params);
if (particle1 != torsionIndexArray[i][0] || particle2 != torsionIndexArray[i][1] || particle3 != torsionIndexArray[i][2] || particle4 != torsionIndexArray[i][3])
throw OpenMMException("updateParametersInContext: The set of particles in a torsion has changed");
for (int j = 0; j < numParameters; j++)
torsionParamArray[i][j] = (RealOpenMM) params[j];
}
}
ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
disposeRealArray(particleParamArray, numParticles);
disposeIntArray(exclusionArray, numParticles);
......
platforms/reference/src/ReferenceKernels.h
View file @ 4bc723ab
......@@ -261,6 +261,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 HarmonicBondForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const HarmonicBondForce& force);
private:
int numBonds;
int **bondIndexArray;
......@@ -291,6 +298,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 CustomBondForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const CustomBondForce& force);
private:
int numBonds;
int **bondIndexArray;
......@@ -323,6 +337,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 HarmonicAngleForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force);
private:
int numAngles;
int **angleIndexArray;
......@@ -353,6 +374,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 CustomAngleForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const CustomAngleForce& force);
private:
int numAngles;
int **angleIndexArray;
......@@ -385,6 +413,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 PeriodicTorsionForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force);
private:
int numTorsions;
int **torsionIndexArray;
......@@ -415,6 +450,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 RBTorsionForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const RBTorsionForce& force);
private:
int numTorsions;
int **torsionIndexArray;
......@@ -474,6 +516,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 CustomTorsionForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const CustomTorsionForce& force);
private:
int numTorsions;
int **torsionIndexArray;
......
platforms/reference/tests/TestReferenceCustomAngleForce.cpp
View file @ 4bc723ab
......@@ -90,11 +90,36 @@ void testAngles() {
init_gen_rand(0, sfmt);
vector<Vec3> positions(4);
for (int i = 0; i < 10; i++) {
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context c1(customSystem, integrator1, platform);
Context c2(harmonicSystem, integrator2, platform);
for (int i = 0; i < 10; i++) {
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = s1.getForces();
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
// Try changing the angle parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[1] = 0.9;
custom->setAngleParameters(0, 0, 1, 2, parameters);
parameters[0] = 2.1;
parameters[1] = 0.6;
custom->setAngleParameters(1, 1, 2, 3, parameters);
custom->updateParametersInContext(c1);
harmonic->setAngleParameters(0, 0, 1, 2, 1.6, 0.9);
harmonic->setAngleParameters(1, 1, 2, 3, 2.1, 0.6);
harmonic->updateParametersInContext(c2);
{
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
......
platforms/reference/tests/TestReferenceCustomBondForce.cpp
View file @ 4bc723ab
......@@ -74,11 +74,31 @@ void testBonds() {
positions[2] = Vec3(1, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, -0.8*0.5, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.7*0.2, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
ASSERT_EQUAL_TOL(0.5*0.8*0.5*0.5 + 0.5*0.7*0.2*0.2, state.getPotentialEnergy(), TOL);
}
// Try changing the bond parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[1] = 0.9;
forceField->setBondParameters(0, 0, 1, parameters);
parameters[0] = 1.3;
parameters[1] = 0.8;
forceField->setBondParameters(1, 1, 2, parameters);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, -0.9*0.4, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.8*0.3, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
ASSERT_EQUAL_TOL(0.5*0.9*0.4*0.4 + 0.5*0.8*0.3*0.3, state.getPotentialEnergy(), TOL);
}
}
int main() {
......
platforms/reference/tests/TestReferenceCustomTorsionForce.cpp
View file @ 4bc723ab
......@@ -111,6 +111,31 @@ void testTorsions() {
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
// Try changing the torsion parameters and make sure it's still correct.
parameters[0] = 1.6;
parameters[1] = 2;
custom->setTorsionParameters(0, 0, 1, 2, 3, parameters);
parameters[0] = 2.1;
parameters[1] = 3;
custom->setTorsionParameters(1, 1, 2, 3, 4, parameters);
custom->updateParametersInContext(c1);
periodic->setTorsionParameters(0, 0, 1, 2, 3, 2, 1.6, 0.5);
periodic->setTorsionParameters(1, 1, 2, 3, 4, 3, 2.1, 0.5);
periodic->updateParametersInContext(c2);
{
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = s1.getForces();
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
}
void testRange() {
......
platforms/reference/tests/TestReferenceHarmonicAngleForce.cpp
View file @ 4bc723ab
......@@ -68,6 +68,7 @@ void testAngles() {
positions[3] = Vec3(2, 1, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double torque1 = 1.1*PI_M/6;
double torque2 = 1.2*PI_M/4;
......@@ -75,6 +76,25 @@ void testAngles() {
ASSERT_EQUAL_VEC(Vec3(-0.5*torque2, 0.5*torque2, 0), forces[3], TOL); // reduced by sqrt(2) due to the bond length, another sqrt(2) due to the angle
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(0.5*1.1*(PI_M/6)*(PI_M/6) + 0.5*1.2*(PI_M/4)*(PI_M/4), state.getPotentialEnergy(), TOL);
}
// Try changing the angle parameters and make sure it's still correct.
forceField->setAngleParameters(0, 0, 1, 2, PI_M/3.1, 1.3);
forceField->setAngleParameters(1, 1, 2, 3, PI_M/2.1, 1.4);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double dtheta1 = (PI_M/2)-(PI_M/3.1);
double dtheta2 = (3*PI_M/4)-(PI_M/2.1);
double torque1 = 1.3*dtheta1;
double torque2 = 1.4*dtheta2;
ASSERT_EQUAL_VEC(Vec3(torque1, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-0.5*torque2, 0.5*torque2, 0), forces[3], TOL); // reduced by sqrt(2) due to the bond length, another sqrt(2) due to the angle
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(0.5*1.3*dtheta1*dtheta1 + 0.5*1.4*dtheta2*dtheta2, state.getPotentialEnergy(), TOL);
}
}
int main() {
......
platforms/reference/tests/TestReferenceHarmonicBondForce.cpp
View file @ 4bc723ab
......@@ -66,11 +66,27 @@ void testBonds() {
positions[2] = Vec3(1, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, -0.8*0.5, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.7*0.2, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
ASSERT_EQUAL_TOL(0.5*0.8*0.5*0.5 + 0.5*0.7*0.2*0.2, state.getPotentialEnergy(), TOL);
}
// Try changing the bond parameters and make sure it's still correct.
forceField->setBondParameters(0, 0, 1, 1.6, 0.9);
forceField->setBondParameters(1, 1, 2, 1.3, 0.8);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0, -0.9*0.4, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.8*0.3, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
ASSERT_EQUAL_TOL(0.5*0.9*0.4*0.4 + 0.5*0.8*0.3*0.3, state.getPotentialEnergy(), TOL);
}
}
int main() {
......
platforms/reference/tests/TestReferencePeriodicTorsionForce.cpp
View file @ 4bc723ab
......@@ -67,12 +67,29 @@ void testPeriodicTorsions() {
positions[3] = Vec3(1, 0, 2);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double torque = -2*1.1*std::sin(2*PI_M/3);
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), TOL);
}
// Try changing the torsion parameters and make sure it's still correct.
forceField->setTorsionParameters(0, 0, 1, 2, 3, 3, PI_M/3.2, 1.3);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double dtheta = (3*PI_M/2)-(PI_M/3.2);
double torque = -3*1.3*std::sin(dtheta);
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
ASSERT_EQUAL_TOL(1.3*(1+std::cos(dtheta)), state.getPotentialEnergy(), TOL);
}
}
int main() {
......
platforms/reference/tests/TestReferenceRBTorsionForce.cpp
View file @ 4bc723ab
......@@ -67,6 +67,7 @@ void testRBTorsions() {
positions[3] = Vec3(1, 1, 1);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double psi = 0.25*PI_M - PI_M;
double torque = 0.0;
......@@ -83,6 +84,31 @@ void testRBTorsions() {
energy += c*std::pow(std::cos(psi), i);
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
// Try changing the torsion parameters and make sure it's still correct.
forceField->setTorsionParameters(0, 0, 1, 2, 3, 0.11, 0.22, 0.33, 0.44, 0.55, 0.66);
forceField->updateParametersInContext(context);
state = context.getState(State::Forces | State::Energy);
{
const vector<Vec3>& forces = state.getForces();
double psi = 0.25*PI_M - PI_M;
double torque = 0.0;
for (int i = 1; i < 6; ++i) {
double c = 0.11*(i+1);
torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
}
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, -0.5*torque), forces[3], TOL);
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
double energy = 0.0;
for (int i = 0; i < 6; ++i) {
double c = 0.11*(i+1);
energy += c*std::pow(std::cos(psi), i);
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
}
int main() {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment