Skip to content

GitLab

Commit c44c956d authored by peastman's avatar peastman
Browse files

Deleted lots of debugging code

parent 41cd79a5
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 150 additions and 555 deletions
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaGeneralizedKirkwoodForce.cpp
View file @ c44c956d
...@@ -275,7 +275,7 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce ...@@ -275,7 +275,7 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
system.addForce(amoebaGeneralizedKirkwoodForce); system.addForce(amoebaGeneralizedKirkwoodForce);
} }
   
static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& forces, double& energy, FILE* log) { static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& forces, double& energy) {
std::vector<Vec3> positions(context.getSystem().getNumParticles()); std::vector<Vec3> positions(context.getSystem().getNumParticles());
   
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03); positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03);
...@@ -296,7 +296,7 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& ...@@ -296,7 +296,7 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>&
// setup for villin // setup for villin
   
static void setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::PolarizationType polarizationType, static void setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::PolarizationType polarizationType,
int includeCavityTerm, std::vector<Vec3>& forces, double& energy, FILE* log) { int includeCavityTerm, std::vector<Vec3>& forces, double& energy) {
   
// beginning of Multipole setup // beginning of Multipole setup
   
...@@ -6980,44 +6980,7 @@ static void setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Polariz ...@@ -6980,44 +6980,7 @@ static void setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Polariz
   
static void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy, static void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy,
const std::vector<Vec3>& expectedForces, const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log) { const std::vector<Vec3>& forces, double tolerance) {
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
if (log) {
double conversion = 1.0/4.184;
double energyAbsDiff = fabs(expectedEnergy - energy);
double energyRelDiff = 2.0*energyAbsDiff/(fabs(expectedEnergy) + fabs(energy) + 1.0e-08);
(void) fprintf(log, "%s: expected energy=%14.7e %14.7e absDiff=%15.7e relDiff=%15.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy,
conversion*energyAbsDiff, conversion*energyRelDiff);
if (conversion != 1.0)conversion *= -0.1;
for (unsigned int ii = 0; ii < forces.size(); ii++) {
double expectedNorm = sqrt(expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2]);
double norm = sqrt(forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2]);
double absDiff = fabs(norm - expectedNorm);
double relDiff = 2.0*absDiff/(fabs(norm) + fabs(expectedNorm) + 1.0e-08);
(void) fprintf(log, "%6u %15.7e %15.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*absDiff, conversion*relDiff,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2], conversion*expectedNorm, conversion*norm);
}
(void) fflush(log);
conversion = 1.0;
(void) fprintf(log, "\n%s: expected energy=%14.7e %14.7e no conversion\n", testName.c_str(), conversion*expectedEnergy, conversion*energy);
if (conversion != 1.0)conversion = -1.0;
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,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2]);
}
(void) fflush(log);
}
#endif
   
for (unsigned int ii = 0; ii < forces.size(); ii++) { for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC_MOD(expectedForces[ii], forces[ii], tolerance, testName); ASSERT_EQUAL_VEC_MOD(expectedForces[ii], forces[ii], tolerance, testName);
...@@ -7029,47 +6992,7 @@ static void compareForcesEnergy(std::string& testName, double expectedEnergy, do ...@@ -7029,47 +6992,7 @@ static void compareForcesEnergy(std::string& testName, double expectedEnergy, do
   
static void compareForceNormsEnergy(std::string& testName, double expectedEnergy, double energy, static void compareForceNormsEnergy(std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces, std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log) { const std::vector<Vec3>& forces, double tolerance) {
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
if (log) {
double conversion = 1.0/4.184;
double energyAbsDiff = fabs(expectedEnergy - energy);
double energyRelDiff = 2.0*energyAbsDiff/(fabs(expectedEnergy) + fabs(energy) + 1.0e-08);
(void) fprintf(log, "%s: expected energy=%14.7e %14.7e absDiff=%15.7e relDiff=%15.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy,
conversion*energyAbsDiff, conversion*energyRelDiff);
if (conversion != 1.0)conversion *= -0.1;
for (unsigned int ii = 0; ii < forces.size(); ii++) {
double expectedNorm = sqrt(expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2]);
double norm = sqrt(forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2]);
double absDiff = fabs((norm - expectedNorm));
double relDiff = 2.0*absDiff/(fabs(norm) + fabs(expectedNorm) + 1.0e-08);
(void) fprintf(log, "%6u %15.7e %15.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e] %15.7e %15.7e\n", ii,
fabs(conversion)*absDiff, relDiff,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2],
fabs(conversion)*expectedNorm, fabs(conversion)*norm);
}
(void) fflush(log);
conversion = 1.0;
(void) fprintf(log, "\n%s: expected energy=%14.7e %14.7e no conversion\n", testName.c_str(), conversion*expectedEnergy, conversion*energy);
if (conversion != 1.0)conversion = -1.0;
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,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2]);
}
(void) fflush(log);
}
#endif
for (unsigned int ii = 0; ii < forces.size(); ii++) { for (unsigned int ii = 0; ii < forces.size(); ii++) {
double expectedNorm = sqrt(expectedForces[ii][0]*expectedForces[ii][0] + double expectedNorm = sqrt(expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] + expectedForces[ii][1]*expectedForces[ii][1] +
...@@ -7098,7 +7021,7 @@ static void compareForceNormsEnergy(std::string& testName, double expectedEnergy ...@@ -7098,7 +7021,7 @@ static void compareForceNormsEnergy(std::string& testName, double expectedEnergy
   
// test GK direct polarization for system comprised of two ammonia molecules // test GK direct polarization for system comprised of two ammonia molecules
   
static void testGeneralizedKirkwoodAmmoniaDirectPolarization(FILE* log) { static void testGeneralizedKirkwoodAmmoniaDirectPolarization() {
   
std::string testName = "testGeneralizedKirkwoodAmmoniaDirectPolarization"; std::string testName = "testGeneralizedKirkwoodAmmoniaDirectPolarization";
   
...@@ -7111,7 +7034,7 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization(FILE* log) { ...@@ -7111,7 +7034,7 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization(FILE* log) {
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0); setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference")); Context context(system, integrator, Platform::getPlatformByName("Reference"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log); getForcesEnergyMultipoleAmmonia(context, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
   
double expectedEnergy = -7.6636680e+01; double expectedEnergy = -7.6636680e+01;
...@@ -7126,12 +7049,12 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization(FILE* log) { ...@@ -7126,12 +7049,12 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization(FILE* log) {
expectedForces[7] = Vec3( 9.7274235e+00, -6.3130458e+01, 1.4949024e+02); expectedForces[7] = Vec3( 9.7274235e+00, -6.3130458e+01, 1.4949024e+02);
   
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
   
// test GK mutual polarization for system comprised of two ammonia molecules // test GK mutual polarization for system comprised of two ammonia molecules
   
static void testGeneralizedKirkwoodAmmoniaMutualPolarization(FILE* log) { static void testGeneralizedKirkwoodAmmoniaMutualPolarization() {
   
std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarization"; std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarization";
   
...@@ -7144,7 +7067,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization(FILE* log) { ...@@ -7144,7 +7067,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization(FILE* log) {
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 0); setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference")); Context context(system, integrator, Platform::getPlatformByName("Reference"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log); getForcesEnergyMultipoleAmmonia(context, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
   
double expectedEnergy = -7.8018875e+01; double expectedEnergy = -7.8018875e+01;
...@@ -7159,13 +7082,13 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization(FILE* log) { ...@@ -7159,13 +7082,13 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization(FILE* log) {
expectedForces[7] = Vec3( 5.3895456e+00, -7.7131137e+01, 1.5826273e+02); expectedForces[7] = Vec3( 5.3895456e+00, -7.7131137e+01, 1.5826273e+02);
   
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
   
// test GK mutual polarization for system comprised of two ammonia molecules // test GK mutual polarization for system comprised of two ammonia molecules
// including cavity term // including cavity term
   
static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(FILE* log) { static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm() {
   
std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm"; std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm";
   
...@@ -7180,7 +7103,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(FILE* ...@@ -7180,7 +7103,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(FILE*
ASSERT(!system.usesPeriodicBoundaryConditions()); ASSERT(!system.usesPeriodicBoundaryConditions());
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference")); Context context(system, integrator, Platform::getPlatformByName("Reference"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log); getForcesEnergyMultipoleAmmonia(context, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
   
double expectedEnergy = -6.0434582e+01; double expectedEnergy = -6.0434582e+01;
...@@ -7195,7 +7118,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(FILE* ...@@ -7195,7 +7118,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(FILE*
expectedForces[7] = Vec3( 7.9205382e+00, -7.3716473e+01, 1.5960993e+02); expectedForces[7] = Vec3( 7.9205382e+00, -7.3716473e+01, 1.5960993e+02);
   
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
// Try changing the particle parameters and make sure it's still correct. // Try changing the particle parameters and make sure it's still correct.
...@@ -7212,7 +7135,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(FILE* ...@@ -7212,7 +7135,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(FILE*
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log); compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance);
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
...@@ -7220,12 +7143,12 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(FILE* ...@@ -7220,12 +7143,12 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(FILE*
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaGeneralizedKirkwoodForce->updateParametersInContext(context); amoebaGeneralizedKirkwoodForce->updateParametersInContext(context);
state1 = context.getState(State::Forces | State::Energy); state1 = context.getState(State::Forces | State::Energy);
compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance, log); compareForcesEnergy(testName, state2.getPotentialEnergy(), state1.getPotentialEnergy(), state2.getForces(), state1.getForces(), tolerance);
} }
   
// test GK direct polarization for villin system // test GK direct polarization for villin system
   
static void testGeneralizedKirkwoodVillinDirectPolarization(FILE* log) { static void testGeneralizedKirkwoodVillinDirectPolarization() {
   
std::string testName = "testGeneralizedKirkwoodVillinDirectPolarization"; std::string testName = "testGeneralizedKirkwoodVillinDirectPolarization";
   
...@@ -7233,7 +7156,7 @@ static void testGeneralizedKirkwoodVillinDirectPolarization(FILE* log) { ...@@ -7233,7 +7156,7 @@ static void testGeneralizedKirkwoodVillinDirectPolarization(FILE* log) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
   
setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Direct, 0, forces, energy, log); setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Direct, 0, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
   
double expectedEnergy = -8.4281157e+03; double expectedEnergy = -8.4281157e+03;
...@@ -7844,12 +7767,12 @@ static void testGeneralizedKirkwoodVillinDirectPolarization(FILE* log) { ...@@ -7844,12 +7767,12 @@ static void testGeneralizedKirkwoodVillinDirectPolarization(FILE* log) {
} }
   
double tolerance = 1.0e-05; double tolerance = 1.0e-05;
compareForceNormsEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForceNormsEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
   
// test GK mutual polarization for villin system // test GK mutual polarization for villin system
   
static void testGeneralizedKirkwoodVillinMutualPolarization(FILE* log) { static void testGeneralizedKirkwoodVillinMutualPolarization() {
   
std::string testName = "testGeneralizedKirkwoodVillinMutualPolarization"; std::string testName = "testGeneralizedKirkwoodVillinMutualPolarization";
   
...@@ -7857,7 +7780,7 @@ static void testGeneralizedKirkwoodVillinMutualPolarization(FILE* log) { ...@@ -7857,7 +7780,7 @@ static void testGeneralizedKirkwoodVillinMutualPolarization(FILE* log) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
   
setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Mutual, 0, forces, energy, log); setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Mutual, 0, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
   
double expectedEnergy = -8.6477811e+03; double expectedEnergy = -8.6477811e+03;
...@@ -8468,7 +8391,7 @@ static void testGeneralizedKirkwoodVillinMutualPolarization(FILE* log) { ...@@ -8468,7 +8391,7 @@ static void testGeneralizedKirkwoodVillinMutualPolarization(FILE* log) {
} }
   
double tolerance = 1.0e-05; double tolerance = 1.0e-05;
compareForceNormsEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForceNormsEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
   
int main(int numberOfArguments, char* argv[]) { int main(int numberOfArguments, char* argv[]) {
...@@ -8477,16 +8400,14 @@ int main(int numberOfArguments, char* argv[]) { ...@@ -8477,16 +8400,14 @@ int main(int numberOfArguments, char* argv[]) {
std::cout << "TestReferenceAmoebaGeneralizedKirkwoodForce running test..." << std::endl; std::cout << "TestReferenceAmoebaGeneralizedKirkwoodForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories(); registerAmoebaReferenceKernelFactories();
   
FILE* log = NULL;
// test direct and mutual polarization cases and // test direct and mutual polarization cases and
// mutual polarization w/ the cavity term // mutual polarization w/ the cavity term
   
testGeneralizedKirkwoodAmmoniaMutualPolarization(log); testGeneralizedKirkwoodAmmoniaMutualPolarization();
testGeneralizedKirkwoodAmmoniaDirectPolarization(log); testGeneralizedKirkwoodAmmoniaDirectPolarization();
testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(log); testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm();
testGeneralizedKirkwoodVillinDirectPolarization(log); testGeneralizedKirkwoodVillinDirectPolarization();
testGeneralizedKirkwoodVillinMutualPolarization(log); testGeneralizedKirkwoodVillinMutualPolarization();
   
} }
catch(const std::exception& e) { catch(const std::exception& e) {
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaInPlaneAngleForce.cpp
View file @ c44c956d
...@@ -79,7 +79,7 @@ static double dotVector3(double* vectorX, double* vectorY) { ...@@ -79,7 +79,7 @@ static double dotVector3(double* vectorX, double* vectorY) {
static void getPrefactorsGivenInPlaneAngleCosine(double cosine, double idealInPlaneAngle, double quadraticK, double cubicK, static void getPrefactorsGivenInPlaneAngleCosine(double cosine, double idealInPlaneAngle, double quadraticK, double cubicK,
double quarticK, double penticK, double sexticK, double quarticK, double penticK, double sexticK,
double* dEdR, double* energyTerm, FILE* log) { double* dEdR, double* energyTerm) {
double angle; double angle;
if (cosine >= 1.0) { if (cosine >= 1.0) {
...@@ -91,12 +91,6 @@ static void getPrefactorsGivenInPlaneAngleCosine(double cosine, double idealInPl ...@@ -91,12 +91,6 @@ static void getPrefactorsGivenInPlaneAngleCosine(double cosine, double idealInPl
else { else {
angle = RADIAN*acos(cosine); angle = RADIAN*acos(cosine);
} }
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "getPrefactorsGivenInPlaneAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealInPlaneAngle);
(void) fflush(log);
}
#endif
double deltaIdeal = angle - idealInPlaneAngle; double deltaIdeal = angle - idealInPlaneAngle;
double deltaIdeal2 = deltaIdeal*deltaIdeal; double deltaIdeal2 = deltaIdeal*deltaIdeal;
...@@ -124,7 +118,7 @@ static void getPrefactorsGivenInPlaneAngleCosine(double cosine, double idealInPl ...@@ -124,7 +118,7 @@ static void getPrefactorsGivenInPlaneAngleCosine(double cosine, double idealInPl
} }
static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& positions, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce, static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& positions, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce,
std::vector<Vec3>& forces, double* energy, FILE* log) { std::vector<Vec3>& forces, double* energy) {
int particle1, particle2, particle3, particle4; int particle1, particle2, particle3, particle4;
double idealInPlaneAngle; double idealInPlaneAngle;
...@@ -135,13 +129,6 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -135,13 +129,6 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
double quarticK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleQuartic(); double quarticK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleQuartic();
double penticK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAnglePentic(); double penticK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAnglePentic();
double sexticK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleSextic(); double sexticK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleSextic();
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaInPlaneAngleForce: bond %d [%d %d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK, cubicK, quarticK, penticK, sexticK);
(void) fflush(log);
}
#endif
// T = AD x CD // T = AD x CD
// P = B + T*delta // P = B + T*delta
...@@ -182,12 +169,6 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -182,12 +169,6 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
double rAp2 = dotVector3(deltaR[AP], deltaR[AP]); double rAp2 = dotVector3(deltaR[AP], deltaR[AP]);
double rCp2 = dotVector3(deltaR[CP], deltaR[CP]); double rCp2 = dotVector3(deltaR[CP], deltaR[CP]);
if (rAp2 <= 0.0 && rCp2 <= 0.0) { if (rAp2 <= 0.0 && rCp2 <= 0.0) {
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaInPlaneAngleForce: rAp2 or rCp2 <= 0.0\n");
(void) fflush(log);
}
#endif
return; return;
} }
...@@ -205,7 +186,7 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -205,7 +186,7 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
double dEdR; double dEdR;
double energyTerm; double energyTerm;
getPrefactorsGivenInPlaneAngleCosine(cosine, idealInPlaneAngle, quadraticK, cubicK, getPrefactorsGivenInPlaneAngleCosine(cosine, idealInPlaneAngle, quadraticK, cubicK,
quarticK, penticK, sexticK, &dEdR, &energyTerm, log); quarticK, penticK, sexticK, &dEdR, &energyTerm);
double termA = -dEdR/(rAp2*rm); double termA = -dEdR/(rAp2*rm);
double termC = dEdR/(rCp2*rm); double termC = dEdR/(rCp2*rm);
...@@ -281,7 +262,7 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -281,7 +262,7 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
} }
static void computeAmoebaInPlaneAngleForces(Context& context, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce, static void computeAmoebaInPlaneAngleForces(Context& context, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
...@@ -297,40 +278,19 @@ static void computeAmoebaInPlaneAngleForces(Context& context, AmoebaInPlaneAngle ...@@ -297,40 +278,19 @@ static void computeAmoebaInPlaneAngleForces(Context& context, AmoebaInPlaneAngle
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for (int ii = 0; ii < AmoebaInPlaneAngleForce.getNumAngles(); ii++) { for (int ii = 0; ii < AmoebaInPlaneAngleForce.getNumAngles(); ii++) {
computeAmoebaInPlaneAngleForce(ii, positions, AmoebaInPlaneAngleForce, expectedForces, expectedEnergy, log); computeAmoebaInPlaneAngleForce(ii, positions, AmoebaInPlaneAngleForce, expectedForces, expectedEnergy);
} }
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaInPlaneAngleForces: expected energy=%14.7e\n", *expectedEnergy);
for (unsigned int ii = 0; ii < positions.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2]);
}
(void) fflush(log);
}
#endif
return;
} }
void compareWithExpectedForceAndEnergy(Context& context, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce, void compareWithExpectedForceAndEnergy(Context& context, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce,
double tolerance, const std::string& idString, FILE* log) { double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaInPlaneAngleForces(context, AmoebaInPlaneAngleForce, expectedForces, &expectedEnergy, log); computeAmoebaInPlaneAngleForces(context, AmoebaInPlaneAngleForce, expectedForces, &expectedEnergy);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaInPlaneAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy());
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++) { for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance); ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
...@@ -338,7 +298,7 @@ void compareWithExpectedForceAndEnergy(Context& context, AmoebaInPlaneAngleForce ...@@ -338,7 +298,7 @@ void compareWithExpectedForceAndEnergy(Context& context, AmoebaInPlaneAngleForce
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance); ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
} }
void testOneAngle(FILE* log) { void testOneAngle() {
System system; System system;
int numberOfParticles = 4; int numberOfParticles = 4;
...@@ -376,7 +336,7 @@ void testOneAngle(FILE* log) { ...@@ -376,7 +336,7 @@ void testOneAngle(FILE* log) {
positions[3] = Vec3(1, 1, 1); positions[3] = Vec3(1, 1, 1);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log); compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
// Try changing the angle parameters and make sure it's still correct. // Try changing the angle parameters and make sure it's still correct.
...@@ -384,14 +344,14 @@ void testOneAngle(FILE* log) { ...@@ -384,14 +344,14 @@ void testOneAngle(FILE* log) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log); compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaInPlaneAngleForce->updateParametersInContext(context); amoebaInPlaneAngleForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log); compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
} }
int main(int numberOfArguments, char* argv[]) { int main(int numberOfArguments, char* argv[]) {
...@@ -399,16 +359,7 @@ int main(int numberOfArguments, char* argv[]) { ...@@ -399,16 +359,7 @@ int main(int numberOfArguments, char* argv[]) {
try { try {
std::cout << "TestReferenceAmoebaInPlaneAngleForce running test..." << std::endl; std::cout << "TestReferenceAmoebaInPlaneAngleForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories(); registerAmoebaReferenceKernelFactories();
FILE* log = NULL; testOneAngle();
//FILE* log = stderr;
//FILE* log = fopen("AmoebaInPlaneAngleForce.log", "w");;
testOneAngle(NULL);
#ifdef AMOEBA_DEBUG
if (log && log != stderr)
(void) fclose(log);
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaOutOfPlaneBendForce.cpp
View file @ c44c956d
...@@ -79,7 +79,7 @@ static double dotVector3(double* vectorX, double* vectorY) { ...@@ -79,7 +79,7 @@ static double dotVector3(double* vectorX, double* vectorY) {
static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce, static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
std::vector<Vec3>& forces, double* energy, FILE* log) { std::vector<Vec3>& forces, double* energy) {
double kAngleCubic = amoebaOutOfPlaneBendForce.getAmoebaGlobalOutOfPlaneBendCubic(); double kAngleCubic = amoebaOutOfPlaneBendForce.getAmoebaGlobalOutOfPlaneBendCubic();
...@@ -90,14 +90,6 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -90,14 +90,6 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
int particle1, particle2, particle3, particle4; int particle1, particle2, particle3, particle4;
double kAngleQuadratic; double kAngleQuadratic;
amoebaOutOfPlaneBendForce.getOutOfPlaneBendParameters(bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic); amoebaOutOfPlaneBendForce.getOutOfPlaneBendParameters(bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic);
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaOutOfPlaneBendForce: bond %d [%d %d %d %d] k=[%10.3e %10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic, kAngleCubic, kAngleQuartic, kAnglePentic, kAngleSextic);
(void) fflush(log);
}
#endif
enum { A, B, C, D, LastAtomIndex }; enum { A, B, C, D, LastAtomIndex };
enum { AB, CB, DB, AD, CD, LastDeltaIndex }; enum { AB, CB, DB, AD, CD, LastDeltaIndex };
...@@ -141,13 +133,6 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -141,13 +133,6 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
else { else {
angle = RADIAN*acos(cosine); angle = RADIAN*acos(cosine);
} }
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaOutOfPlaneBendForce: bkk2=%14.7e rDB2=%14.7e cos=%14.7e dt=%14.7e]\n",
bkk2, rDB2, cosine, angle);
(void) fflush(log);
}
#endif
// chain rule // chain rule
...@@ -244,7 +229,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -244,7 +229,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
} }
static void computeAmoebaOutOfPlaneBendForces(Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce, static void computeAmoebaOutOfPlaneBendForces(Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
...@@ -260,41 +245,19 @@ static void computeAmoebaOutOfPlaneBendForces(Context& context, AmoebaOutOfPlane ...@@ -260,41 +245,19 @@ static void computeAmoebaOutOfPlaneBendForces(Context& context, AmoebaOutOfPlane
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for (int ii = 0; ii < amoebaOutOfPlaneBendForce.getNumOutOfPlaneBends(); ii++) { for (int ii = 0; ii < amoebaOutOfPlaneBendForce.getNumOutOfPlaneBends(); ii++) {
computeAmoebaOutOfPlaneBendForce(ii, positions, amoebaOutOfPlaneBendForce, expectedForces, expectedEnergy, log); computeAmoebaOutOfPlaneBendForce(ii, positions, amoebaOutOfPlaneBendForce, expectedForces, expectedEnergy);
}
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e\n", *expectedEnergy);
for (unsigned int ii = 0; ii < positions.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2]);
}
(void) fflush(log);
} }
#endif
return;
} }
void compareWithExpectedForceAndEnergy(Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce, void compareWithExpectedForceAndEnergy(Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
double tolerance, const std::string& idString, FILE* log) { double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaOutOfPlaneBendForces(context, amoebaOutOfPlaneBendForce, expectedForces, &expectedEnergy, log); computeAmoebaOutOfPlaneBendForces(context, amoebaOutOfPlaneBendForce, expectedForces, &expectedEnergy);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy());
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++) { for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance); ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
...@@ -302,7 +265,7 @@ void compareWithExpectedForceAndEnergy(Context& context, AmoebaOutOfPlaneBendFor ...@@ -302,7 +265,7 @@ void compareWithExpectedForceAndEnergy(Context& context, AmoebaOutOfPlaneBendFor
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance); ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
} }
void testOneOutOfPlaneBend(FILE* log) { void testOneOutOfPlaneBend() {
System system; System system;
int numberOfParticles = 4; int numberOfParticles = 4;
...@@ -336,7 +299,7 @@ void testOneOutOfPlaneBend(FILE* log) { ...@@ -336,7 +299,7 @@ void testOneOutOfPlaneBend(FILE* log) {
positions[3] = Vec3(0.245568230E+02, 0.250215290E+02, 0.796852800E+01); positions[3] = Vec3(0.245568230E+02, 0.250215290E+02, 0.796852800E+01);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log); compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
// Try changing the bend parameters and make sure it's still correct. // Try changing the bend parameters and make sure it's still correct.
...@@ -344,17 +307,17 @@ void testOneOutOfPlaneBend(FILE* log) { ...@@ -344,17 +307,17 @@ void testOneOutOfPlaneBend(FILE* log) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log); compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaOutOfPlaneBendForce->updateParametersInContext(context); amoebaOutOfPlaneBendForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log); compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
} }
void testOneOutOfPlaneBend2(FILE* log, int setId) { void testOneOutOfPlaneBend2(int setId) {
System system; System system;
int numberOfParticles = 4; int numberOfParticles = 4;
...@@ -444,20 +407,10 @@ void testOneOutOfPlaneBend2(FILE* log, int setId) { ...@@ -444,20 +407,10 @@ void testOneOutOfPlaneBend2(FILE* log, int setId) {
kOutOfPlaneBend = 0.214755281E-01; kOutOfPlaneBend = 0.214755281E-01;
} }
else { else {
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "Set id %d not recognized.\n", setId);
}
#endif
std::stringstream buffer; std::stringstream buffer;
buffer << "Set id " << setId << " not recognized."; buffer << "Set id " << setId << " not recognized.";
throw OpenMMException(buffer.str()); throw OpenMMException(buffer.str());
} }
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "Set id %d.\n", setId);
}
#endif
amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend); amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend);
...@@ -467,11 +420,6 @@ void testOneOutOfPlaneBend2(FILE* log, int setId) { ...@@ -467,11 +420,6 @@ void testOneOutOfPlaneBend2(FILE* log, int setId) {
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numberOfParticles); ii++) { for (unsigned int ii = 0; ii < static_cast<unsigned int>(numberOfParticles); ii++) {
if (coordinates.find(particleIndices[ii]) == coordinates.end()) { if (coordinates.find(particleIndices[ii]) == coordinates.end()) {
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "Coordinates %d not loaded.", particleIndices[ii]);
}
#endif
std::stringstream buffer; std::stringstream buffer;
buffer << "Coordinates " << particleIndices[ii] << " not loaded."; buffer << "Coordinates " << particleIndices[ii] << " not loaded.";
throw OpenMMException(buffer.str()); throw OpenMMException(buffer.str());
...@@ -480,7 +428,7 @@ void testOneOutOfPlaneBend2(FILE* log, int setId) { ...@@ -480,7 +428,7 @@ void testOneOutOfPlaneBend2(FILE* log, int setId) {
} }
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log); compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
static int iter = 0; static int iter = 0;
static std::map<int,Vec3> totalForces; static std::map<int,Vec3> totalForces;
...@@ -503,7 +451,7 @@ void testOneOutOfPlaneBend2(FILE* log, int setId) { ...@@ -503,7 +451,7 @@ void testOneOutOfPlaneBend2(FILE* log, int setId) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
forces.resize(numberOfParticles); forces.resize(numberOfParticles);
double energy; double energy;
computeAmoebaOutOfPlaneBendForce(0, positions, *amoebaOutOfPlaneBendForce, forces, &energy, log); computeAmoebaOutOfPlaneBendForce(0, positions, *amoebaOutOfPlaneBendForce, forces, &energy);
totalEnergy += energy; totalEnergy += energy;
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numberOfParticles); ii++) { for (unsigned int ii = 0; ii < static_cast<unsigned int>(numberOfParticles); ii++) {
...@@ -511,22 +459,6 @@ void testOneOutOfPlaneBend2(FILE* log, int setId) { ...@@ -511,22 +459,6 @@ void testOneOutOfPlaneBend2(FILE* log, int setId) {
totalForces[particleIndices[ii]][jj] += forces[ii][jj]; totalForces[particleIndices[ii]][jj] += forces[ii][jj];
} }
} }
if (iter == 6) {
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaOutOfPlaneBendForces: energy=%14.7e\n", totalEnergy);
for (std::map<int,Vec3>::iterator ii = totalForces.begin(); ii != totalForces.end(); ii++) {
int particleIndex = ii->first;
Vec3 forces = ii->second;
(void) fprintf(log, "%6d [%14.7e %14.7e %14.7e] \n", particleIndex,
forces[0], forces[1], forces[2]);
}
(void) fflush(log);
}
#endif
}
} }
int main(int numberOfArguments, char* argv[]) { int main(int numberOfArguments, char* argv[]) {
...@@ -534,21 +466,11 @@ int main(int numberOfArguments, char* argv[]) { ...@@ -534,21 +466,11 @@ int main(int numberOfArguments, char* argv[]) {
try { try {
std::cout << "TestReferenceAmoebaOutOfPlaneBendForce running test..." << std::endl; std::cout << "TestReferenceAmoebaOutOfPlaneBendForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories(); registerAmoebaReferenceKernelFactories();
testOneOutOfPlaneBend();
//FILE* log = stderr; //testOneOutOfPlaneBend2(atoi(argv[1]));
FILE* log = NULL;
//FILE* log = fopen("AmoebaOutOfPlaneBendForce.log", "w");;
testOneOutOfPlaneBend(log);
//testOneOutOfPlaneBend2(log, atoi(argv[1]));
//for (int ii = 1; ii <= 6; ii++) { //for (int ii = 1; ii <= 6; ii++) {
// testOneOutOfPlaneBend2(log, ii); // testOneOutOfPlaneBend2(ii);
//} //}
#ifdef AMOEBA_DEBUG
if (log && log != stderr)
(void) fclose(log);
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaPiTorsionForce.cpp
View file @ c44c956d
...@@ -78,19 +78,12 @@ static double dotVector3(double* vectorX, double* vectorY) { ...@@ -78,19 +78,12 @@ static double dotVector3(double* vectorX, double* vectorY) {
static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& positions, AmoebaPiTorsionForce& amoebaPiTorsionForce, static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& positions, AmoebaPiTorsionForce& amoebaPiTorsionForce,
std::vector<Vec3>& forces, double* energy, FILE* log) { std::vector<Vec3>& forces, double* energy) {
int particle1, particle2, particle3, particle4, particle5, particle6; int particle1, particle2, particle3, particle4, particle5, particle6;
double kTorsion; double kTorsion;
amoebaPiTorsionForce.getPiTorsionParameters(bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion); amoebaPiTorsionForce.getPiTorsionParameters(bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion);
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaPiTorsionForce: bond %d [%d %d %d %d %d %d] k=%10.3e\n",
bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion);
(void) fflush(log);
}
#endif
enum { AD, BD, EC, FC, P, Q, CP, DC, QD, T, U, TU, DP, QC, dT, dU, dP, dQ, dC1, dC2, dD1, dD2, LastDeltaIndex }; enum { AD, BD, EC, FC, P, Q, CP, DC, QD, T, U, TU, DP, QC, dT, dU, dP, dQ, dC1, dC2, dD1, dD2, LastDeltaIndex };
double deltaR[LastDeltaIndex][3]; double deltaR[LastDeltaIndex][3];
...@@ -212,7 +205,7 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit ...@@ -212,7 +205,7 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
} }
static void computeAmoebaPiTorsionForces(Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce, static void computeAmoebaPiTorsionForces(Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
...@@ -228,40 +221,19 @@ static void computeAmoebaPiTorsionForces(Context& context, AmoebaPiTorsionForce& ...@@ -228,40 +221,19 @@ static void computeAmoebaPiTorsionForces(Context& context, AmoebaPiTorsionForce&
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for (int ii = 0; ii < amoebaPiTorsionForce.getNumPiTorsions(); ii++) { for (int ii = 0; ii < amoebaPiTorsionForce.getNumPiTorsions(); ii++) {
computeAmoebaPiTorsionForce(ii, positions, amoebaPiTorsionForce, expectedForces, expectedEnergy, log); computeAmoebaPiTorsionForce(ii, positions, amoebaPiTorsionForce, expectedForces, expectedEnergy);
}
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaPiTorsionForces: expected energy=%14.7e\n", *expectedEnergy);
for (unsigned int ii = 0; ii < positions.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2]);
}
(void) fflush(log);
} }
#endif
return;
} }
void compareWithExpectedForceAndEnergy(Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce, void compareWithExpectedForceAndEnergy(Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
double tolerance, const std::string& idString, FILE* log) { double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaPiTorsionForces(context, amoebaPiTorsionForce, expectedForces, &expectedEnergy, log); computeAmoebaPiTorsionForces(context, amoebaPiTorsionForce, expectedForces, &expectedEnergy);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaPiTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy());
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++) { for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance); ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
...@@ -269,7 +241,7 @@ void compareWithExpectedForceAndEnergy(Context& context, AmoebaPiTorsionForce& a ...@@ -269,7 +241,7 @@ void compareWithExpectedForceAndEnergy(Context& context, AmoebaPiTorsionForce& a
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance); ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
} }
void testOnePiTorsion(FILE* log) { void testOnePiTorsion() {
System system; System system;
int numberOfParticles = 6; int numberOfParticles = 6;
...@@ -300,7 +272,7 @@ void testOnePiTorsion(FILE* log) { ...@@ -300,7 +272,7 @@ void testOnePiTorsion(FILE* log) {
positions[5] = Vec3(0.254124630E+02, 0.234691880E+02, 0.773335400E+01); positions[5] = Vec3(0.254124630E+02, 0.234691880E+02, 0.773335400E+01);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log); compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion");
// Try changing the torsion parameters and make sure it's still correct. // Try changing the torsion parameters and make sure it's still correct.
...@@ -308,14 +280,14 @@ void testOnePiTorsion(FILE* log) { ...@@ -308,14 +280,14 @@ void testOnePiTorsion(FILE* log) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log); compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaPiTorsionForce->updateParametersInContext(context); amoebaPiTorsionForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log); compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion");
} }
int main(int numberOfArguments, char* argv[]) { int main(int numberOfArguments, char* argv[]) {
...@@ -323,16 +295,7 @@ int main(int numberOfArguments, char* argv[]) { ...@@ -323,16 +295,7 @@ int main(int numberOfArguments, char* argv[]) {
try { try {
std::cout << "TestReferenceAmoebaPiTorsionForce running test..." << std::endl; std::cout << "TestReferenceAmoebaPiTorsionForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories(); registerAmoebaReferenceKernelFactories();
FILE* log = NULL; testOnePiTorsion();
//FILE* log = stderr;
//FILE* log = fopen("AmoebaPiTorsionForce1.log", "w");;
testOnePiTorsion(log);
#ifdef AMOEBA_DEBUG
if (log && log != stderr)
(void) fclose(log);
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaStretchBendForce.cpp
View file @ c44c956d
...@@ -80,20 +80,13 @@ static double dotVector3(double* vectorX, double* vectorY) { ...@@ -80,20 +80,13 @@ static double dotVector3(double* vectorX, double* vectorY) {
static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaStretchBendForce& amoebaStretchBendForce, static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaStretchBendForce& amoebaStretchBendForce,
std::vector<Vec3>& forces, double* energy, FILE* log) { std::vector<Vec3>& forces, double* energy) {
int particle1, particle2, particle3; int particle1, particle2, particle3;
double abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend; double abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend;
amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend); amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend);
angleStretchBend *= RADIAN; angleStretchBend *= RADIAN;
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaStretchBendForce: bond %d [%d %d %d] ab=%10.3e cb=%10.3e angle=%10.3e k1=%10.3e k2=%10.3e\n",
bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend);
(void) fflush(log);
}
#endif
enum { A, B, C, LastAtomIndex }; enum { A, B, C, LastAtomIndex };
enum { AB, CB, CBxAB, ABxP, CBxP, LastDeltaIndex }; enum { AB, CB, CBxAB, ABxP, CBxP, LastDeltaIndex };
...@@ -189,18 +182,10 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -189,18 +182,10 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
forces[particle3][2] -= subForce[2][2]; forces[particle3][2] -= subForce[2][2];
*energy += dt*drkk; *energy += dt*drkk;
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaStretchBendForce: angle=%10.3e dt=%10.3e dr=%10.3e\n", angle, dt, dr);
(void) fflush(log);
}
#endif
return;
} }
static void computeAmoebaStretchBendForces(Context& context, AmoebaStretchBendForce& amoebaStretchBendForce, static void computeAmoebaStretchBendForces(Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
...@@ -216,48 +201,26 @@ static void computeAmoebaStretchBendForces(Context& context, AmoebaStretchBendFo ...@@ -216,48 +201,26 @@ static void computeAmoebaStretchBendForces(Context& context, AmoebaStretchBendFo
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for (int ii = 0; ii < amoebaStretchBendForce.getNumStretchBends(); ii++) { for (int ii = 0; ii < amoebaStretchBendForce.getNumStretchBends(); ii++) {
computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy, log); computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy);
}
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaStretchBendForces: expected energy=%14.7e\n", *expectedEnergy);
for (unsigned int ii = 0; ii < positions.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2]);
}
(void) fflush(log);
} }
#endif
return;
} }
void compareWithExpectedForceAndEnergy(Context& context, AmoebaStretchBendForce& amoebaStretchBendForce, void compareWithExpectedForceAndEnergy(Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
double tolerance, const std::string& idString, FILE* log) { double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaStretchBendForces(context, amoebaStretchBendForce, expectedForces, &expectedEnergy, log); computeAmoebaStretchBendForces(context, amoebaStretchBendForce, expectedForces, &expectedEnergy);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaStretchBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy());
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++) { for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance); ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
} }
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance); ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
} }
void testOneStretchBend(FILE* log) { void testOneStretchBend() {
System system; System system;
int numberOfParticles = 3; int numberOfParticles = 3;
...@@ -289,7 +252,7 @@ void testOneStretchBend(FILE* log) { ...@@ -289,7 +252,7 @@ void testOneStretchBend(FILE* log) {
positions[2] = Vec3(0.269573220E+02, 0.236108860E+02, 0.216376800E+01); positions[2] = Vec3(0.269573220E+02, 0.236108860E+02, 0.216376800E+01);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log); compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
// Try changing the stretch-bend parameters and make sure it's still correct. // Try changing the stretch-bend parameters and make sure it's still correct.
...@@ -297,14 +260,14 @@ void testOneStretchBend(FILE* log) { ...@@ -297,14 +260,14 @@ void testOneStretchBend(FILE* log) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log); compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaStretchBendForce->updateParametersInContext(context); amoebaStretchBendForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log); compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
} }
int main(int numberOfArguments, char* argv[]) { int main(int numberOfArguments, char* argv[]) {
...@@ -312,16 +275,7 @@ int main(int numberOfArguments, char* argv[]) { ...@@ -312,16 +275,7 @@ int main(int numberOfArguments, char* argv[]) {
try { try {
std::cout << "TestReferenceAmoebaStretchBendForce running test..." << std::endl; std::cout << "TestReferenceAmoebaStretchBendForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories(); registerAmoebaReferenceKernelFactories();
testOneStretchBend();
FILE* log = NULL;
//FILE* log = stderr;
//FILE* log = fopen("AmoebaStretchBendForce1.log", "w");;
testOneStretchBend(log);
#ifdef AMOEBA_DEBUG
if (log && log != stderr)
(void) fclose(log);
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaVdwForce.cpp
View file @ c44c956d
...@@ -57,7 +57,7 @@ extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories(); ...@@ -57,7 +57,7 @@ extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();
const double TOL = 1e-4; const double TOL = 1e-4;
void testVdw(FILE* log) { void testVdw() {
System system; System system;
int numberOfParticles = 6; int numberOfParticles = 6;
...@@ -162,17 +162,6 @@ void testVdw(FILE* log) { ...@@ -162,17 +162,6 @@ void testVdw(FILE* log) {
} }
expectedEnergy *= CalToJoule; expectedEnergy *= CalToJoule;
#ifdef AMOEBA_DEBUG
if (log) {
(void) fprintf(log, "computeAmoebaVdwForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy());
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
double tolerance = 1.0e-03; double tolerance = 1.0e-03;
for (unsigned int ii = 0; ii < forces.size(); ii++) { for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance); ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
...@@ -210,7 +199,7 @@ void testVdw(FILE* log) { ...@@ -210,7 +199,7 @@ void testVdw(FILE* log) {
} }
void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff, void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff,
double boxDimension, std::vector<Vec3>& forces, double& energy, FILE* log) { double boxDimension, std::vector<Vec3>& forces, double& energy) {
// beginning of Vdw setup // beginning of Vdw setup
...@@ -353,24 +342,7 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co ...@@ -353,24 +342,7 @@ void setupAndGetForcesEnergyVdwAmmonia(const std::string& sigmaCombiningRule, co
void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy, void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces, std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance, FILE* log) { std::vector<Vec3>& forces, double tolerance) {
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
if (log) {
double conversion = 1.0/4.184;
(void) fprintf(log, "%s: expected energy=%14.7e %14.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy);
conversion *= -0.1;
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,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2]);
}
(void) fflush(log);
}
#endif
for (unsigned int ii = 0; ii < forces.size(); ii++) { for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC_MOD(expectedForces[ii], forces[ii], tolerance, testName); ASSERT_EQUAL_VEC_MOD(expectedForces[ii], forces[ii], tolerance, testName);
} }
...@@ -379,7 +351,7 @@ void compareForcesEnergy(std::string& testName, double expectedEnergy, double en ...@@ -379,7 +351,7 @@ void compareForcesEnergy(std::string& testName, double expectedEnergy, double en
// test VDW w/ sigmaRule=CubicMean and epsilonRule=HHG // test VDW w/ sigmaRule=CubicMean and epsilonRule=HHG
void testVdwAmmoniaCubicMeanHhg(FILE* log) { void testVdwAmmoniaCubicMeanHhg() {
std::string testName = "testVdwAmmoniaCubicMeanHhg"; std::string testName = "testVdwAmmoniaCubicMeanHhg";
...@@ -389,7 +361,7 @@ void testVdwAmmoniaCubicMeanHhg(FILE* log) { ...@@ -389,7 +361,7 @@ void testVdwAmmoniaCubicMeanHhg(FILE* log) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia("CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy, log); setupAndGetForcesEnergyVdwAmmonia("CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 4.8012258e+00; double expectedEnergy = 4.8012258e+00;
...@@ -404,12 +376,12 @@ void testVdwAmmoniaCubicMeanHhg(FILE* log) { ...@@ -404,12 +376,12 @@ void testVdwAmmoniaCubicMeanHhg(FILE* log) {
expectedForces[7] = Vec3( 1.6756544e+00, 3.2497316e-01, -1.7906832e-01); expectedForces[7] = Vec3( 1.6756544e+00, 3.2497316e-01, -1.7906832e-01);
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
// test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic // test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
void testVdwAmmoniaArithmeticArithmetic(FILE* log) { void testVdwAmmoniaArithmeticArithmetic() {
std::string testName = "testVdwAmmoniaArithmeticArithmetic"; std::string testName = "testVdwAmmoniaArithmeticArithmetic";
...@@ -419,7 +391,7 @@ void testVdwAmmoniaArithmeticArithmetic(FILE* log) { ...@@ -419,7 +391,7 @@ void testVdwAmmoniaArithmeticArithmetic(FILE* log) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia("ARITHMETIC", "ARITHMETIC", cutoff, boxDimension, forces, energy, log); setupAndGetForcesEnergyVdwAmmonia("ARITHMETIC", "ARITHMETIC", cutoff, boxDimension, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 4.2252403e+00; double expectedEnergy = 4.2252403e+00;
...@@ -434,12 +406,12 @@ void testVdwAmmoniaArithmeticArithmetic(FILE* log) { ...@@ -434,12 +406,12 @@ void testVdwAmmoniaArithmeticArithmetic(FILE* log) {
expectedForces[7] = Vec3( 2.3761408e+00, 4.6871961e-01, -2.4731607e-01); expectedForces[7] = Vec3( 2.3761408e+00, 4.6871961e-01, -2.4731607e-01);
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
// test VDW w/ sigmaRule=Geometric and epsilonRule=Geometric // test VDW w/ sigmaRule=Geometric and epsilonRule=Geometric
void testVdwAmmoniaGeometricGeometric(FILE* log) { void testVdwAmmoniaGeometricGeometric() {
std::string testName = "testVdwAmmoniaGeometricGeometric"; std::string testName = "testVdwAmmoniaGeometricGeometric";
...@@ -448,7 +420,7 @@ void testVdwAmmoniaGeometricGeometric(FILE* log) { ...@@ -448,7 +420,7 @@ void testVdwAmmoniaGeometricGeometric(FILE* log) {
double cutoff = 9000000.0; double cutoff = 9000000.0;
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia("GEOMETRIC", "GEOMETRIC", cutoff, boxDimension, forces, energy, log); setupAndGetForcesEnergyVdwAmmonia("GEOMETRIC", "GEOMETRIC", cutoff, boxDimension, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 2.5249914e+00; double expectedEnergy = 2.5249914e+00;
...@@ -463,10 +435,10 @@ void testVdwAmmoniaGeometricGeometric(FILE* log) { ...@@ -463,10 +435,10 @@ void testVdwAmmoniaGeometricGeometric(FILE* log) {
expectedForces[7] = Vec3( 1.8109211e+00, 3.5273117e-01, -1.9224723e-01); expectedForces[7] = Vec3( 1.8109211e+00, 3.5273117e-01, -1.9224723e-01);
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
void testVdwAmmoniaCubicMeanHarmonic(FILE* log) { void testVdwAmmoniaCubicMeanHarmonic() {
std::string testName = "testVdwAmmoniaCubicMeanHarmonic"; std::string testName = "testVdwAmmoniaCubicMeanHarmonic";
...@@ -475,7 +447,7 @@ void testVdwAmmoniaCubicMeanHarmonic(FILE* log) { ...@@ -475,7 +447,7 @@ void testVdwAmmoniaCubicMeanHarmonic(FILE* log) {
double cutoff = 9000000.0; double cutoff = 9000000.0;
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia("CUBIC-MEAN", "HARMONIC", cutoff, boxDimension, forces, energy, log); setupAndGetForcesEnergyVdwAmmonia("CUBIC-MEAN", "HARMONIC", cutoff, boxDimension, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 4.1369069e+00; double expectedEnergy = 4.1369069e+00;
...@@ -490,14 +462,14 @@ void testVdwAmmoniaCubicMeanHarmonic(FILE* log) { ...@@ -490,14 +462,14 @@ void testVdwAmmoniaCubicMeanHarmonic(FILE* log) {
expectedForces[7] = Vec3( 1.5080748e+00, 2.9058422e-01, -1.6274118e-01); expectedForces[7] = Vec3( 1.5080748e+00, 2.9058422e-01, -1.6274118e-01);
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
// test w/ cutoff=0.25 nm; single ixn between two particles (0 and 6); force nonzero on // test w/ cutoff=0.25 nm; single ixn between two particles (0 and 6); force nonzero on
// particle 4 due to reduction applied to NH // particle 4 due to reduction applied to NH
// the distance between 0 and 6 is ~ 0.235 so the ixn is in the tapered region // the distance between 0 and 6 is ~ 0.235 so the ixn is in the tapered region
void testVdwTaper(FILE* log) { void testVdwTaper() {
std::string testName = "testVdwTaper"; std::string testName = "testVdwTaper";
...@@ -507,7 +479,7 @@ void testVdwTaper(FILE* log) { ...@@ -507,7 +479,7 @@ void testVdwTaper(FILE* log) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia("CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy, log); setupAndGetForcesEnergyVdwAmmonia("CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 3.5478444e+00; double expectedEnergy = 3.5478444e+00;
...@@ -522,12 +494,12 @@ void testVdwTaper(FILE* log) { ...@@ -522,12 +494,12 @@ void testVdwTaper(FILE* log) {
expectedForces[7] = Vec3( -0.0000000e+00, -0.0000000e+00, -0.0000000e+00); expectedForces[7] = Vec3( -0.0000000e+00, -0.0000000e+00, -0.0000000e+00);
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
// test PBC // test PBC
void testVdwPBC(FILE* log) { void testVdwPBC() {
std::string testName = "testVdwPBC"; std::string testName = "testVdwPBC";
...@@ -537,7 +509,7 @@ void testVdwPBC(FILE* log) { ...@@ -537,7 +509,7 @@ void testVdwPBC(FILE* log) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia("CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy, log); setupAndGetForcesEnergyVdwAmmonia("CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 1.4949141e+01; double expectedEnergy = 1.4949141e+01;
...@@ -555,14 +527,14 @@ void testVdwPBC(FILE* log) { ...@@ -555,14 +527,14 @@ void testVdwPBC(FILE* log) {
// if tapering turned off, then absolute difference < 2.0e-05 // if tapering turned off, then absolute difference < 2.0e-05
double tolerance = 5.0e-04; double tolerance = 5.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
// create box of 216 water molecules // create box of 216 water molecules
void setupAndGetForcesEnergyVdwWater(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff, void setupAndGetForcesEnergyVdwWater(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff,
double boxDimension, int includeVdwDispersionCorrection, double boxDimension, int includeVdwDispersionCorrection,
std::vector<Vec3>& forces, double& energy, FILE* log) { std::vector<Vec3>& forces, double& energy) {
// beginning of Vdw setup // beginning of Vdw setup
...@@ -1282,7 +1254,7 @@ void setupAndGetForcesEnergyVdwWater(const std::string& sigmaCombiningRule, cons ...@@ -1282,7 +1254,7 @@ void setupAndGetForcesEnergyVdwWater(const std::string& sigmaCombiningRule, cons
// test employing box of 216 water molecules w/ and w/o dispersion correction // test employing box of 216 water molecules w/ and w/o dispersion correction
void testVdwWater(int includeVdwDispersionCorrection, FILE* log) { void testVdwWater(int includeVdwDispersionCorrection) {
std::string testName; std::string testName;
...@@ -1299,7 +1271,7 @@ void testVdwWater(int includeVdwDispersionCorrection, FILE* log) { ...@@ -1299,7 +1271,7 @@ void testVdwWater(int includeVdwDispersionCorrection, FILE* log) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwWater("CUBIC-MEAN", "HHG", cutoff, boxDimension, includeVdwDispersionCorrection, forces, energy, log); setupAndGetForcesEnergyVdwWater("CUBIC-MEAN", "HHG", cutoff, boxDimension, includeVdwDispersionCorrection, forces, energy);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
// initialize expected energy and forces // initialize expected energy and forces
...@@ -1965,7 +1937,7 @@ void testVdwWater(int includeVdwDispersionCorrection, FILE* log) { ...@@ -1965,7 +1937,7 @@ void testVdwWater(int includeVdwDispersionCorrection, FILE* log) {
// if tapering turned off, then absolute difference < 2.0e-05 // if tapering turned off, then absolute difference < 2.0e-05
double tolerance = 5.0e-04; double tolerance = 5.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log); compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
// test sigma/epsilon rules for dispersion correction // test sigma/epsilon rules for dispersion correction
...@@ -1988,7 +1960,7 @@ void testVdwWater(int includeVdwDispersionCorrection, FILE* log) { ...@@ -1988,7 +1960,7 @@ void testVdwWater(int includeVdwDispersionCorrection, FILE* log) {
expectedEnergies.push_back(3.2774624e+03); expectedEnergies.push_back(3.2774624e+03);
for (unsigned int ii = 0; ii < sigmaRules.size(); ii++) { for (unsigned int ii = 0; ii < sigmaRules.size(); ii++) {
setupAndGetForcesEnergyVdwWater(sigmaRules[ii], epsilonRules[ii], cutoff, boxDimension, includeVdwDispersionCorrection, forces, energy, log); setupAndGetForcesEnergyVdwWater(sigmaRules[ii], epsilonRules[ii], cutoff, boxDimension, includeVdwDispersionCorrection, forces, energy);
testName = "testVdwWaterWithDispersionCorrection_" + sigmaRules[ii] + '_' + epsilonRules[ii]; testName = "testVdwWaterWithDispersionCorrection_" + sigmaRules[ii] + '_' + epsilonRules[ii];
ASSERT_EQUAL_TOL_MOD(expectedEnergies[ii], energy, tolerance, testName); ASSERT_EQUAL_TOL_MOD(expectedEnergies[ii], energy, tolerance, testName);
} }
...@@ -2057,49 +2029,46 @@ int main(int numberOfArguments, char* argv[]) { ...@@ -2057,49 +2029,46 @@ int main(int numberOfArguments, char* argv[]) {
try { try {
std::cout << "TestReferenceAmoebaVdwForce running test..." << std::endl; std::cout << "TestReferenceAmoebaVdwForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories(); registerAmoebaReferenceKernelFactories();
testVdw();
FILE* log = NULL;
testVdw(log);
// tests using two ammonia molecules // tests using two ammonia molecules
// test VDW w/ sigmaRule=CubicMean and epsilonRule=HHG // test VDW w/ sigmaRule=CubicMean and epsilonRule=HHG
testVdwAmmoniaCubicMeanHhg(log); testVdwAmmoniaCubicMeanHhg();
// test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic // test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
testVdwAmmoniaArithmeticArithmetic(log); testVdwAmmoniaArithmeticArithmetic();
// test VDW w/ sigmaRule=Geometric and epsilonRule=Geometric // test VDW w/ sigmaRule=Geometric and epsilonRule=Geometric
testVdwAmmoniaGeometricGeometric(log); testVdwAmmoniaGeometricGeometric();
// test VDW w/ sigmaRule=CubicMean and epsilonRule=Harmonic // test VDW w/ sigmaRule=CubicMean and epsilonRule=Harmonic
testVdwAmmoniaCubicMeanHarmonic(log); testVdwAmmoniaCubicMeanHarmonic();
// test w/ cutoff=0.25 nm; single ixn between two particles (0 and 6); force nonzero on // test w/ cutoff=0.25 nm; single ixn between two particles (0 and 6); force nonzero on
// particle 4 due to reduction applied to NH // particle 4 due to reduction applied to NH
// the distance between 0 and 6 is ~ 0.235 so the ixn is in the tapered region // the distance between 0 and 6 is ~ 0.235 so the ixn is in the tapered region
testVdwTaper(log); testVdwTaper();
// test PBC // test PBC
testVdwPBC(log); testVdwPBC();
// tests based on box of water // tests based on box of water
int includeVdwDispersionCorrection = 0; int includeVdwDispersionCorrection = 0;
testVdwWater(includeVdwDispersionCorrection, log); testVdwWater(includeVdwDispersionCorrection);
// includes tests for various combinations of sigma/epsilon rules // includes tests for various combinations of sigma/epsilon rules
// when computing vdw dispersion correction // when computing vdw dispersion correction
includeVdwDispersionCorrection = 1; includeVdwDispersionCorrection = 1;
testVdwWater(includeVdwDispersionCorrection, log); testVdwWater(includeVdwDispersionCorrection);
// test triclinic boxes // test triclinic boxes
......
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