Skip to content

GitLab

Commit 162d69a2 authored by peastman's avatar peastman
Browse files

Merge pull request #827 from peastman/reference 

Lots of cleanup to reference platform
parents 2379b982 c44c956d
Hide whitespace changes
Inline Side-by-side
Showing with 3941 additions and 4375 deletions
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaBondForce.cpp
View file @ 162d69a2
...@@ -49,22 +49,22 @@ extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories(); ...@@ -49,22 +49,22 @@ extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();
const double TOL = 1e-5; const double TOL = 1e-5;
static void computeAmoebaBondForce(int bondIndex, std::vector<Vec3>& positions, AmoebaBondForce& AmoebaBondForce, static void computeAmoebaBondForce(int bondIndex, std::vector<Vec3>& positions, AmoebaBondForce& AmoebaBondForce,
std::vector<Vec3>& forces, double* energy ) { std::vector<Vec3>& forces, double* energy) {
int particle1, particle2; int particle1, particle2;
double bondLength; double bondLength;
double quadraticK; double quadraticK;
double cubicK = AmoebaBondForce.getAmoebaGlobalBondCubic(); double cubicK = AmoebaBondForce.getAmoebaGlobalBondCubic();
double quarticK = AmoebaBondForce.getAmoebaGlobalBondQuartic(); double quarticK = AmoebaBondForce.getAmoebaGlobalBondQuartic();
AmoebaBondForce.getBondParameters(bondIndex, particle1, particle2, bondLength, quadraticK ); AmoebaBondForce.getBondParameters(bondIndex, particle1, particle2, bondLength, quadraticK);
double deltaR[3]; double deltaR[3];
double r2 = 0.0; double r2 = 0.0;
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[ii] = positions[particle2][ii] - positions[particle1][ii]; deltaR[ii] = positions[particle2][ii] - positions[particle1][ii];
r2 += deltaR[ii]*deltaR[ii]; r2 += deltaR[ii]*deltaR[ii];
} }
double r = sqrt( r2 ); double r = sqrt(r2);
double bondDelta = (r - bondLength); double bondDelta = (r - bondLength);
double bondDelta2 = bondDelta*bondDelta; double bondDelta2 = bondDelta*bondDelta;
...@@ -84,64 +84,43 @@ static void computeAmoebaBondForce(int bondIndex, std::vector<Vec3>& positions, ...@@ -84,64 +84,43 @@ static void computeAmoebaBondForce(int bondIndex, std::vector<Vec3>& positions,
} }
static void computeAmoebaBondForces( Context& context, AmoebaBondForce& AmoebaBondForce, static void computeAmoebaBondForces(Context& context, AmoebaBondForce& AmoebaBondForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
State state = context.getState(State::Positions); State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions(); std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() ); expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){ for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0; expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
} }
// calculates forces/energy // calculates forces/energy
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for( int ii = 0; ii < AmoebaBondForce.getNumBonds(); ii++ ){ for (int ii = 0; ii < AmoebaBondForce.getNumBonds(); ii++) {
computeAmoebaBondForce(ii, positions, AmoebaBondForce, expectedForces, expectedEnergy ); computeAmoebaBondForce(ii, positions, AmoebaBondForce, expectedForces, expectedEnergy);
} }
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaBondForces: expected energy=%15.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%15.7e %15.7e %15.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
}
(void) fflush( log );
}
#endif
return;
} }
void compareWithExpectedForceAndEnergy( Context& context, AmoebaBondForce& AmoebaBondForce, double tolerance, const std::string& idString, FILE* log) { void compareWithExpectedForceAndEnergy(Context& context, AmoebaBondForce& AmoebaBondForce, double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaBondForces( context, AmoebaBondForce, expectedForces, &expectedEnergy, NULL ); computeAmoebaBondForces(context, AmoebaBondForce, 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, "computeAmoebaBondForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.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 testOneBond( FILE* log ) { void testOneBond() {
System system; System system;
...@@ -156,22 +135,22 @@ void testOneBond( FILE* log ) { ...@@ -156,22 +135,22 @@ void testOneBond( FILE* log ) {
double quadraticK = 1.0; double quadraticK = 1.0;
double cubicK = 2.0; double cubicK = 2.0;
double quarticicK = 3.0; double quarticicK = 3.0;
amoebaBondForce->setAmoebaGlobalBondCubic( cubicK ); amoebaBondForce->setAmoebaGlobalBondCubic(cubicK);
amoebaBondForce->setAmoebaGlobalBondQuartic( quarticicK ); amoebaBondForce->setAmoebaGlobalBondQuartic(quarticicK);
amoebaBondForce->addBond(0, 1, bondLength, quadraticK); amoebaBondForce->addBond(0, 1, bondLength, quadraticK);
system.addForce(amoebaBondForce); system.addForce(amoebaBondForce);
Context context(system, integrator, Platform::getPlatformByName( "Reference")); Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(2); std::vector<Vec3> positions(2);
positions[0] = Vec3(0, 1, 0); positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0); positions[1] = Vec3(0, 0, 0);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testOneBond", log ); compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testOneBond");
} }
void testTwoBond( FILE* log ) { void testTwoBond() {
System system; System system;
...@@ -187,15 +166,15 @@ void testTwoBond( FILE* log ) { ...@@ -187,15 +166,15 @@ void testTwoBond( FILE* log ) {
double quadraticK = 1.0; double quadraticK = 1.0;
double cubicK = 2.0; double cubicK = 2.0;
double quarticicK = 3.0; double quarticicK = 3.0;
amoebaBondForce->setAmoebaGlobalBondCubic( cubicK ); amoebaBondForce->setAmoebaGlobalBondCubic(cubicK);
amoebaBondForce->setAmoebaGlobalBondQuartic( quarticicK ); amoebaBondForce->setAmoebaGlobalBondQuartic(quarticicK);
amoebaBondForce->addBond(0, 1, bondLength, quadraticK); amoebaBondForce->addBond(0, 1, bondLength, quadraticK);
amoebaBondForce->addBond(1, 2, bondLength, quadraticK); amoebaBondForce->addBond(1, 2, bondLength, quadraticK);
system.addForce(amoebaBondForce); system.addForce(amoebaBondForce);
ASSERT(!amoebaBondForce->usesPeriodicBoundaryConditions()); ASSERT(!amoebaBondForce->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions()); ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, Platform::getPlatformByName( "Reference")); Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(3); std::vector<Vec3> positions(3);
positions[0] = Vec3(0, 1, 0); positions[0] = Vec3(0, 1, 0);
...@@ -203,7 +182,7 @@ void testTwoBond( FILE* log ) { ...@@ -203,7 +182,7 @@ void testTwoBond( FILE* log ) {
positions[2] = Vec3(1, 0, 1); positions[2] = Vec3(1, 0, 1);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log ); compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond");
// Try changing the bond parameters and make sure it's still correct. // Try changing the bond parameters and make sure it's still correct.
...@@ -212,31 +191,23 @@ void testTwoBond( FILE* log ) { ...@@ -212,31 +191,23 @@ void testTwoBond( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log ); compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaBondForce->updateParametersInContext(context); amoebaBondForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log ); compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond");
} }
int main( int numberOfArguments, char* argv[] ) { int main(int numberOfArguments, char* argv[]) {
try { try {
std::cout << "TestReferenceAmoebaBondForce running test..." << std::endl; std::cout << "TestReferenceAmoebaBondForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories(); registerAmoebaReferenceKernelFactories();
FILE* log = NULL; //testOneBond();
//FILE* log = stderr; testTwoBond();
//testOneBond( log );
testTwoBond( 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/TestReferenceAmoebaGeneralizedKirkwoodForce.cpp
View file @ 162d69a2
...@@ -65,10 +65,10 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce ...@@ -65,10 +65,10 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();; AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
int numberOfParticles = 8; int numberOfParticles = 8;
   
amoebaMultipoleForce->setNonbondedMethod( AmoebaMultipoleForce::NoCutoff ); amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::NoCutoff);
amoebaMultipoleForce->setPolarizationType( polarizationType ); amoebaMultipoleForce->setPolarizationType(polarizationType);
amoebaMultipoleForce->setMutualInducedTargetEpsilon( 1.0e-06 ); amoebaMultipoleForce->setMutualInducedTargetEpsilon(1.0e-06);
amoebaMultipoleForce->setMutualInducedMaxIterations( 500 ); amoebaMultipoleForce->setMutualInducedMaxIterations(500);
   
std::vector<double> nitrogenMolecularDipole(3); std::vector<double> nitrogenMolecularDipole(3);
std::vector<double> nitrogenMolecularQuadrupole(9); std::vector<double> nitrogenMolecularQuadrupole(9);
...@@ -89,8 +89,8 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce ...@@ -89,8 +89,8 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
   
// first N // first N
   
system.addParticle( 1.4007000e+01 ); system.addParticle(1.4007000e+01);
amoebaMultipoleForce->addMultipole( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 1, 2, 3, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03 ); amoebaMultipoleForce->addMultipole( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 1, 2, 3, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03);
   
// 3 H attached to first N // 3 H attached to first N
   
...@@ -110,182 +110,182 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce ...@@ -110,182 +110,182 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
hydrogenMolecularQuadrupole[7] = 0.0000000e+00; hydrogenMolecularQuadrupole[7] = 0.0000000e+00;
hydrogenMolecularQuadrupole[8] = 2.4549167e-06; hydrogenMolecularQuadrupole[8] = 2.4549167e-06;
   
system.addParticle( 1.0080000e+00 ); system.addParticle(1.0080000e+00);
system.addParticle( 1.0080000e+00 ); system.addParticle(1.0080000e+00);
system.addParticle( 1.0080000e+00 ); system.addParticle(1.0080000e+00);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 2, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 2, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 2, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 2, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
   
// second N // second N
   
system.addParticle( 1.4007000e+01 ); system.addParticle( 1.4007000e+01);
amoebaMultipoleForce->addMultipole( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 5, 6, 7, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03 ); amoebaMultipoleForce->addMultipole( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 5, 6, 7, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03);
   
// 3 H attached to second N // 3 H attached to second N
   
system.addParticle( 1.0080000e+00 ); system.addParticle( 1.0080000e+00);
system.addParticle( 1.0080000e+00 ); system.addParticle( 1.0080000e+00);
system.addParticle( 1.0080000e+00 ); system.addParticle( 1.0080000e+00);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 6, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 6, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 6, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 ); amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 6, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
   
// covalent maps // covalent maps
   
std::vector< int > covalentMap; std::vector< int > covalentMap;
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 0 ); covalentMap.push_back(0);
covalentMap.push_back( 1 ); covalentMap.push_back(1);
covalentMap.push_back( 2 ); covalentMap.push_back(2);
covalentMap.push_back( 3 ); covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap ); amoebaMultipoleForce->setCovalentMap(7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap ); amoebaMultipoleForce->setCovalentMap(7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
   
covalentMap.resize(0); covalentMap.resize(0);
covalentMap.push_back( 4 ); covalentMap.push_back(4);
covalentMap.push_back( 5 ); covalentMap.push_back(5);
covalentMap.push_back( 6 ); covalentMap.push_back(6);
covalentMap.push_back( 7 ); covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap ); amoebaMultipoleForce->setCovalentMap(7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
   
system.addForce(amoebaMultipoleForce); system.addForce(amoebaMultipoleForce);
   
// GK force // GK force
   
amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01 ); amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01);
amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00 ); amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00);
amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm( includeCavityTerm ); amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm(includeCavityTerm);
   
// addParticle: charge, radius, scalingFactor // addParticle: charge, radius, scalingFactor
   
for( unsigned int ii = 0; ii < 2; ii++ ){ for (unsigned int ii = 0; ii < 2; ii++) {
amoebaGeneralizedKirkwoodForce->addParticle( -5.7960000e-01, 1.5965000e-01, 6.9000000e-01 ); amoebaGeneralizedKirkwoodForce->addParticle( -5.7960000e-01, 1.5965000e-01, 6.9000000e-01);
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 ); amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01);
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 ); amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01);
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 ); amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01);
} }
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);
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02 ); positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02);
positions[2] = Vec3( 2.0843610e-01, 8.0953200e-02, 3.7462200e-02 ); positions[2] = Vec3( 2.0843610e-01, 8.0953200e-02, 3.7462200e-02);
positions[3] = Vec3( 1.7280780e-01, 2.0730000e-04, -9.8741700e-02 ); positions[3] = Vec3( 1.7280780e-01, 2.0730000e-04, -9.8741700e-02);
positions[4] = Vec3( -1.6743680e-01, 1.5900000e-05, -6.6149000e-03 ); positions[4] = Vec3( -1.6743680e-01, 1.5900000e-05, -6.6149000e-03);
positions[5] = Vec3( -2.0428260e-01, 8.1071500e-02, 4.1343900e-02 ); positions[5] = Vec3( -2.0428260e-01, 8.1071500e-02, 4.1343900e-02);
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02 ); positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02);
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02 ); positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02);
   
context.setPositions(positions); context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
...@@ -295,8 +295,8 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& ...@@ -295,8 +295,8 @@ 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
   
...@@ -305,13 +305,13 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -305,13 +305,13 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();; AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
int numberOfParticles = 596; int numberOfParticles = 596;
   
amoebaMultipoleForce->setNonbondedMethod( AmoebaMultipoleForce::NoCutoff ); amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::NoCutoff);
amoebaMultipoleForce->setPolarizationType( polarizationType ); amoebaMultipoleForce->setPolarizationType(polarizationType);
amoebaMultipoleForce->setMutualInducedTargetEpsilon( 1.0e-06 ); amoebaMultipoleForce->setMutualInducedTargetEpsilon(1.0e-06);
amoebaMultipoleForce->setMutualInducedMaxIterations( 500 ); amoebaMultipoleForce->setMutualInducedMaxIterations(500);
   
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle( 1.0 ); system.addParticle(1.0);
} }
   
static const double multipoleData[] = { static const double multipoleData[] = {
...@@ -929,7 +929,7 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -929,7 +929,7 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
std::vector<double> quadrupole(9); std::vector<double> quadrupole(9);
unsigned int entriesPerParticle = 21; unsigned int entriesPerParticle = 21;
const double* data = multipoleData; const double* data = multipoleData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
   
dipole[0] = data[dipoleIndex + 0]; dipole[0] = data[dipoleIndex + 0];
dipole[1] = data[dipoleIndex + 1]; dipole[1] = data[dipoleIndex + 1];
...@@ -945,9 +945,9 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -945,9 +945,9 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
quadrupole[7] = data[quadrupoleIndex + 7]; quadrupole[7] = data[quadrupoleIndex + 7];
quadrupole[8] = data[quadrupoleIndex + 8]; quadrupole[8] = data[quadrupoleIndex + 8];
   
amoebaMultipoleForce->addMultipole( data[chargeIndex], dipole, quadrupole, static_cast<int>(data[axisTypeIndex]), amoebaMultipoleForce->addMultipole(data[chargeIndex], dipole, quadrupole, static_cast<int>(data[axisTypeIndex]),
static_cast<int>(data[multipoleAtomZIndex]), static_cast<int>(data[multipoleAtomXIndex]), static_cast<int>(data[multipoleAtomYIndex]), static_cast<int>(data[multipoleAtomZIndex]), static_cast<int>(data[multipoleAtomXIndex]), static_cast<int>(data[multipoleAtomYIndex]),
data[tholeIndex], data[dampingFactorIndex], data[polarityIndex] ); data[tholeIndex], data[dampingFactorIndex], data[polarityIndex]);
data += entriesPerParticle; data += entriesPerParticle;
} }
   
...@@ -5727,15 +5727,15 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -5727,15 +5727,15 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
unsigned int covalentMapDataSize = sizeof(covalentMapData)/sizeof(int); unsigned int covalentMapDataSize = sizeof(covalentMapData)/sizeof(int);
   
unsigned int covalentIndex = 0; unsigned int covalentIndex = 0;
while( covalentIndex < covalentMapDataSize ){ while (covalentIndex < covalentMapDataSize) {
int particleIndex = covalentMapData[covalentIndex++]; int particleIndex = covalentMapData[covalentIndex++];
int typeIndex = covalentMapData[covalentIndex++]; int typeIndex = covalentMapData[covalentIndex++];
int entries = covalentMapData[covalentIndex++]; int entries = covalentMapData[covalentIndex++];
std::vector< int > covalentMap(entries); std::vector< int > covalentMap(entries);
for( unsigned int ii = 0; ii < entries; ii++ ){ for (unsigned int ii = 0; ii < entries; ii++) {
covalentMap[ii] = covalentMapData[covalentIndex++]; covalentMap[ii] = covalentMapData[covalentIndex++];
} }
amoebaMultipoleForce->setCovalentMap( particleIndex, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(typeIndex), covalentMap ); amoebaMultipoleForce->setCovalentMap(particleIndex, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(typeIndex), covalentMap);
} }
system.addForce(amoebaMultipoleForce); system.addForce(amoebaMultipoleForce);
   
...@@ -5744,9 +5744,9 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -5744,9 +5744,9 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
// GK force // GK force
   
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce(); AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01 ); amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01);
amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00 ); amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00);
amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm( includeCavityTerm ); amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm(includeCavityTerm);
   
// addParticle: charge, radius, scalingFactor // addParticle: charge, radius, scalingFactor
   
...@@ -6350,8 +6350,8 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -6350,8 +6350,8 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
}; };
   
const double* gkData = generalizedKirkwoodData; const double* gkData = generalizedKirkwoodData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
amoebaGeneralizedKirkwoodForce->addParticle( gkData[0], gkData[1], gkData[2] ); amoebaGeneralizedKirkwoodForce->addParticle(gkData[0], gkData[1], gkData[2]);
gkData += 3; gkData += 3;
} }
system.addForce(amoebaGeneralizedKirkwoodForce); system.addForce(amoebaGeneralizedKirkwoodForce);
...@@ -6960,15 +6960,15 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -6960,15 +6960,15 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
}; };
   
const double* positionsDataPtr = positionsData; const double* positionsDataPtr = positionsData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
positions[ii] = Vec3( positionsDataPtr[0], positionsDataPtr[1], positionsDataPtr[2] ); positions[ii] = Vec3(positionsDataPtr[0], positionsDataPtr[1], positionsDataPtr[2]);
positionsDataPtr += 3; positionsDataPtr += 3;
} }
   
std::string platformName; std::string platformName;
platformName = "Reference"; platformName = "Reference";
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName( platformName ) ); Context context(system, integrator, Platform::getPlatformByName(platformName));
   
context.setPositions(positions); context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
...@@ -6978,117 +6978,40 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -6978,117 +6978,40 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
   
// compare forces and energies // compare forces and energies
   
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);
} }
ASSERT_EQUAL_TOL_MOD( expectedEnergy, energy, tolerance, testName ); ASSERT_EQUAL_TOL_MOD(expectedEnergy, energy, tolerance, testName);
} }
   
// compare relative differences in force norms and energies // compare relative differences in force norms and energies
   
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) {
for (unsigned int ii = 0; ii < forces.size(); ii++) {
double expectedNorm = sqrt(expectedForces[ii][0]*expectedForces[ii][0] +
//#define AMOEBA_DEBUG expectedForces[ii][1]*expectedForces[ii][1] +
#ifdef AMOEBA_DEBUG expectedForces[ii][2]*expectedForces[ii][2]);
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++ ){
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 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 absDiff = fabs(norm - expectedNorm);
double relDiff = 2.0*absDiff/(fabs( norm ) + fabs( expectedNorm ) + 1.0e-08); double relDiff = 2.0*absDiff/(fabs(norm) + fabs(expectedNorm) + 1.0e-08);
   
if( relDiff > tolerance && absDiff > 0.001 ){ if (relDiff > tolerance && absDiff > 0.001) {
std::stringstream details; std::stringstream details;
details << testName << "Relative difference in norms " << relDiff << " larger than allowed tolerance at particle=" << ii; details << testName << "Relative difference in norms " << relDiff << " larger than allowed tolerance at particle=" << ii;
details << ": norms=" << norm << " expected norm=" << expectedNorm; details << ": norms=" << norm << " expected norm=" << expectedNorm;
throwException(__FILE__, __LINE__, details.str()); throwException(__FILE__, __LINE__, details.str());
} }
} }
double energyAbsDiff = fabs( expectedEnergy - energy ); double energyAbsDiff = fabs(expectedEnergy - energy);
double energyRelDiff = 2.0*energyAbsDiff/( fabs( expectedEnergy ) + fabs( energy ) + 1.0e-08 ); double energyRelDiff = 2.0*energyAbsDiff/(fabs(expectedEnergy) + fabs(energy) + 1.0e-08);
if( energyRelDiff > tolerance ){ if (energyRelDiff > tolerance) {
std::stringstream details; std::stringstream details;
details << testName << "Relative difference in energies " << energyRelDiff << " larger than allowed tolerance."; details << testName << "Relative difference in energies " << energyRelDiff << " larger than allowed tolerance.";
details << "Energies=" << energy << " expected energy=" << expectedEnergy; details << "Energies=" << energy << " expected energy=" << expectedEnergy;
...@@ -7098,7 +7021,7 @@ static void compareForceNormsEnergy( std::string& testName, double expectedEnerg ...@@ -7098,7 +7021,7 @@ static void compareForceNormsEnergy( std::string& testName, double expectedEnerg
   
// 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,27 +7034,27 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) { ...@@ -7111,27 +7034,27 @@ 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;
   
expectedForces[0] = Vec3( -6.9252994e+02, -8.9085133e+00, 9.6489739e+01 ); expectedForces[0] = Vec3( -6.9252994e+02, -8.9085133e+00, 9.6489739e+01);
expectedForces[1] = Vec3( 1.5593797e+02, -6.0331931e+01, 1.5104507e+01 ); expectedForces[1] = Vec3( 1.5593797e+02, -6.0331931e+01, 1.5104507e+01);
expectedForces[2] = Vec3( 1.5870088e+02, 6.1702809e+01, 6.7708985e+00 ); expectedForces[2] = Vec3( 1.5870088e+02, 6.1702809e+01, 6.7708985e+00);
expectedForces[3] = Vec3( 1.4089885e+02, 7.5870617e+00, -1.1362294e+02 ); expectedForces[3] = Vec3( 1.4089885e+02, 7.5870617e+00, -1.1362294e+02);
expectedForces[4] = Vec3( -1.8916205e+02, 2.1465549e-01, -4.3433152e+02 ); expectedForces[4] = Vec3( -1.8916205e+02, 2.1465549e-01, -4.3433152e+02);
expectedForces[5] = Vec3( 1.0208290e+01, 6.2676753e+01, 1.4987953e+02 ); expectedForces[5] = Vec3( 1.0208290e+01, 6.2676753e+01, 1.4987953e+02);
expectedForces[6] = Vec3( 4.0621859e+02, 1.8962203e-01, 1.3021956e+02 ); expectedForces[6] = Vec3( 4.0621859e+02, 1.8962203e-01, 1.3021956e+02);
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,28 +7067,28 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) { ...@@ -7144,28 +7067,28 @@ 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;
   
expectedForces[0] = Vec3( -7.6820301e+02, -1.0102760e+01, 1.0094389e+02 ); expectedForces[0] = Vec3( -7.6820301e+02, -1.0102760e+01, 1.0094389e+02);
expectedForces[1] = Vec3( 1.7037307e+02, -7.5621857e+01, 2.3320365e+01 ); expectedForces[1] = Vec3( 1.7037307e+02, -7.5621857e+01, 2.3320365e+01);
expectedForces[2] = Vec3( 1.7353828e+02, 7.7199741e+01, 1.3965379e+01 ); expectedForces[2] = Vec3( 1.7353828e+02, 7.7199741e+01, 1.3965379e+01);
expectedForces[3] = Vec3( 1.5045244e+02, 8.5784569e+00, -1.3377619e+02 ); expectedForces[3] = Vec3( 1.5045244e+02, 8.5784569e+00, -1.3377619e+02);
expectedForces[4] = Vec3( -2.1811615e+02, -1.6818022e-01, -4.6103163e+02 ); expectedForces[4] = Vec3( -2.1811615e+02, -1.6818022e-01, -4.6103163e+02);
expectedForces[5] = Vec3( 6.2091942e+00, 7.6748687e+01, 1.5883463e+02 ); expectedForces[5] = Vec3( 6.2091942e+00, 7.6748687e+01, 1.5883463e+02);
expectedForces[6] = Vec3( 4.8035662e+02, 4.9704902e-01, 1.3948083e+02 ); expectedForces[6] = Vec3( 4.8035662e+02, 4.9704902e-01, 1.3948083e+02);
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,22 +7103,22 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE ...@@ -7180,22 +7103,22 @@ 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;
   
expectedForces[0] = Vec3( -7.8323218e+02, -1.0097644e+01, 1.0256890e+02 ); expectedForces[0] = Vec3( -7.8323218e+02, -1.0097644e+01, 1.0256890e+02);
expectedForces[1] = Vec3( 1.7078480e+02, -7.1896701e+01, 2.0840172e+01 ); expectedForces[1] = Vec3( 1.7078480e+02, -7.1896701e+01, 2.0840172e+01);
expectedForces[2] = Vec3( 1.7394089e+02, 7.3488594e+01, 1.1484648e+01 ); expectedForces[2] = Vec3( 1.7394089e+02, 7.3488594e+01, 1.1484648e+01);
expectedForces[3] = Vec3( 1.5169364e+02, 8.5611299e+00, -1.2968050e+02 ); expectedForces[3] = Vec3( 1.5169364e+02, 8.5611299e+00, -1.2968050e+02);
expectedForces[4] = Vec3( -2.1669693e+02, -1.5926823e-01, -4.6636274e+02 ); expectedForces[4] = Vec3( -2.1669693e+02, -1.5926823e-01, -4.6636274e+02);
expectedForces[5] = Vec3( 8.7397444e+00, 7.3330990e+01, 1.6016898e+02 ); expectedForces[5] = Vec3( 8.7397444e+00, 7.3330990e+01, 1.6016898e+02);
expectedForces[6] = Vec3( 4.8684950e+02, 4.8937161e-01, 1.4137061e+02 ); expectedForces[6] = Vec3( 4.8684950e+02, 4.8937161e-01, 1.4137061e+02);
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;
...@@ -7838,18 +7761,18 @@ static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) { ...@@ -7838,18 +7761,18 @@ static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) {
}; };
   
const double* forceDataPtr = forceData; const double* forceDataPtr = forceData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
expectedForces[ii] = Vec3( forceDataPtr[0], forceDataPtr[1], forceDataPtr[2] ); expectedForces[ii] = Vec3(forceDataPtr[0], forceDataPtr[1], forceDataPtr[2]);
forceDataPtr += 3; forceDataPtr += 3;
} }
   
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;
...@@ -8462,31 +8385,29 @@ static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) { ...@@ -8462,31 +8385,29 @@ static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) {
}; };
   
const double* forceDataPtr = forceData; const double* forceDataPtr = forceData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
expectedForces[ii] = Vec3( forceDataPtr[0], forceDataPtr[1], forceDataPtr[2] ); expectedForces[ii] = Vec3(forceDataPtr[0], forceDataPtr[1], forceDataPtr[2]);
forceDataPtr += 3; forceDataPtr += 3;
} }
   
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[]) {
   
try { try {
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 @ 162d69a2
...@@ -64,7 +64,7 @@ const double TOL = 1e-5; ...@@ -64,7 +64,7 @@ const double TOL = 1e-5;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){ static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1]; vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2]; vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
...@@ -73,30 +73,24 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto ...@@ -73,30 +73,24 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return; return;
} }
static double dotVector3( double* vectorX, double* vectorY ){ static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2]; return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
} }
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) {
angle = 0.0f; angle = 0.0f;
} }
else if( cosine <= -1.0 ){ else if (cosine <= -1.0) {
angle = RADIAN*PI_M; angle = RADIAN*PI_M;
} }
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;
...@@ -105,11 +99,11 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP ...@@ -105,11 +99,11 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP
// deltaIdeal = r - r_0 // deltaIdeal = r - r_0
*dEdR = ( 2.0 + *dEdR = (2.0 +
3.0*cubicK* deltaIdeal + 3.0*cubicK* deltaIdeal +
4.0*quarticK*deltaIdeal2 + 4.0*quarticK*deltaIdeal2 +
5.0*penticK* deltaIdeal3 + 5.0*penticK* deltaIdeal3 +
6.0*sexticK* deltaIdeal4 ); 6.0*sexticK* deltaIdeal4 );
*dEdR *= RADIAN*quadraticK*deltaIdeal; *dEdR *= RADIAN*quadraticK*deltaIdeal;
...@@ -124,24 +118,17 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP ...@@ -124,24 +118,17 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP
} }
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;
double quadraticK; double quadraticK;
AmoebaInPlaneAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK ); AmoebaInPlaneAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK);
double cubicK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleCubic(); double cubicK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleCubic();
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
...@@ -162,82 +149,76 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -162,82 +149,76 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
// APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex // APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex
double deltaR[LastDeltaAtomIndex][3]; double deltaR[LastDeltaAtomIndex][3];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii]; deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii];
deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii]; deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii];
deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii]; deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii];
} }
crossProductVector3( deltaR[AD], deltaR[CD], deltaR[T] ); crossProductVector3(deltaR[AD], deltaR[CD], deltaR[T]);
double rT2 = dotVector3( deltaR[T], deltaR[T] ); double rT2 = dotVector3(deltaR[T], deltaR[T]);
double delta = dotVector3( deltaR[T], deltaR[BD] ); double delta = dotVector3(deltaR[T], deltaR[BD]);
delta *= -1.0/rT2; delta *= -1.0/rT2;
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[P][ii] = positions[particle2][ii] + deltaR[T][ii]*delta; deltaR[P][ii] = positions[particle2][ii] + deltaR[T][ii]*delta;
deltaR[AP][ii] = positions[particle1][ii] - deltaR[P][ii]; deltaR[AP][ii] = positions[particle1][ii] - deltaR[P][ii];
deltaR[CP][ii] = positions[particle3][ii] - deltaR[P][ii]; deltaR[CP][ii] = positions[particle3][ii] - deltaR[P][ii];
} }
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;
} }
crossProductVector3( deltaR[CP], deltaR[AP], deltaR[M] ); crossProductVector3(deltaR[CP], deltaR[AP], deltaR[M]);
double rm = dotVector3( deltaR[M], deltaR[M] ); double rm = dotVector3(deltaR[M], deltaR[M]);
rm = sqrt( rm ); rm = sqrt(rm);
if( rm < 0.000001 ){ if (rm < 0.000001) {
rm = 0.000001; rm = 0.000001;
} }
double dot = dotVector3( deltaR[AP], deltaR[CP] ); double dot = dotVector3(deltaR[AP], deltaR[CP]);
double cosine = dot/sqrt( rAp2*rCp2 ); double cosine = dot/sqrt(rAp2*rCp2);
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);
crossProductVector3( deltaR[AP], deltaR[M], deltaR[APxM] ); crossProductVector3(deltaR[AP], deltaR[M], deltaR[APxM]);
crossProductVector3( deltaR[CP], deltaR[M], deltaR[CPxM] ); crossProductVector3(deltaR[CP], deltaR[M], deltaR[CPxM]);
// forces will be gathered here // forces will be gathered here
enum { dA, dB, dC, dD, LastDIndex }; enum { dA, dB, dC, dD, LastDIndex };
double forceTerm[LastDIndex][3]; double forceTerm[LastDIndex][3];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
forceTerm[dA][ii] = deltaR[APxM][ii]*termA; forceTerm[dA][ii] = deltaR[APxM][ii]*termA;
forceTerm[dC][ii] = deltaR[CPxM][ii]*termC; forceTerm[dC][ii] = deltaR[CPxM][ii]*termC;
forceTerm[dB][ii] = -1.0*( forceTerm[dA][ii] + forceTerm[dC][ii] ); forceTerm[dB][ii] = -1.0*(forceTerm[dA][ii] + forceTerm[dC][ii]);
} }
double pTrT2 = dotVector3( forceTerm[dB], deltaR[T] ); double pTrT2 = dotVector3(forceTerm[dB], deltaR[T]);
pTrT2 /= rT2; pTrT2 /= rT2;
crossProductVector3( deltaR[CD], forceTerm[dB], deltaR[CDxdB] ); crossProductVector3(deltaR[CD], forceTerm[dB], deltaR[CDxdB]);
crossProductVector3( forceTerm[dB], deltaR[AD], deltaR[dBxAD] ); crossProductVector3(forceTerm[dB], deltaR[AD], deltaR[dBxAD]);
if( fabs( pTrT2 ) > 1.0e-08 ){ if (fabs(pTrT2) > 1.0e-08) {
double delta2 = delta*2.0; double delta2 = delta*2.0;
crossProductVector3( deltaR[BD], deltaR[CD], deltaR[BDxCD] ); crossProductVector3(deltaR[BD], deltaR[CD], deltaR[BDxCD]);
crossProductVector3( deltaR[T], deltaR[CD], deltaR[TxCD] ); crossProductVector3(deltaR[T], deltaR[CD], deltaR[TxCD] );
crossProductVector3( deltaR[AD], deltaR[BD], deltaR[ADxBD] ); crossProductVector3(deltaR[AD], deltaR[BD], deltaR[ADxBD]);
crossProductVector3( deltaR[AD], deltaR[T], deltaR[ADxT] ); crossProductVector3(deltaR[AD], deltaR[T], deltaR[ADxT] );
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
double term = deltaR[BDxCD][ii] + delta2*deltaR[TxCD][ii]; double term = deltaR[BDxCD][ii] + delta2*deltaR[TxCD][ii];
forceTerm[dA][ii] += delta*deltaR[CDxdB][ii] + term*pTrT2; forceTerm[dA][ii] += delta*deltaR[CDxdB][ii] + term*pTrT2;
...@@ -245,16 +226,16 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -245,16 +226,16 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
term = deltaR[ADxBD][ii] + delta2*deltaR[ADxT][ii]; term = deltaR[ADxBD][ii] + delta2*deltaR[ADxT][ii];
forceTerm[dC][ii] += delta*deltaR[dBxAD][ii] + term*pTrT2; forceTerm[dC][ii] += delta*deltaR[dBxAD][ii] + term*pTrT2;
forceTerm[dD][ii] = -( forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii] ); forceTerm[dD][ii] = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
} }
} }
else { else {
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
forceTerm[dA][ii] += delta*deltaR[CDxdB][ii]; forceTerm[dA][ii] += delta*deltaR[CDxdB][ii];
forceTerm[dC][ii] += delta*deltaR[dBxAD][ii]; forceTerm[dC][ii] += delta*deltaR[dBxAD][ii];
forceTerm[dD][ii] = -( forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii] ); forceTerm[dD][ii] = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
} }
} }
...@@ -280,69 +261,48 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -280,69 +261,48 @@ 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
State state = context.getState(State::Positions); State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions(); std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() ); expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){ for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0; expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
} }
// calculates forces/energy // calculates forces/energy
*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);
} }
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;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
...@@ -366,7 +326,7 @@ void testOneAngle( FILE* log ) { ...@@ -366,7 +326,7 @@ void testOneAngle( FILE* log ) {
system.addForce(amoebaInPlaneAngleForce); system.addForce(amoebaInPlaneAngleForce);
ASSERT(!amoebaInPlaneAngleForce->usesPeriodicBoundaryConditions()); ASSERT(!amoebaInPlaneAngleForce->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions()); ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, Platform::getPlatformByName( "Reference")); Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
...@@ -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,31 +344,22 @@ void testOneAngle( FILE* log ) { ...@@ -384,31 +344,22 @@ 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[]) {
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/TestReferenceAmoebaMultipoleForce.cpp
View file @ 162d69a2
This source diff could not be displayed because it is too large. You can view the blob instead.
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaOutOfPlaneBendForce.cpp
View file @ 162d69a2
...@@ -64,7 +64,7 @@ const double TOL = 1e-3; ...@@ -64,7 +64,7 @@ const double TOL = 1e-3;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){ static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1]; vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2]; vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
...@@ -73,13 +73,13 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto ...@@ -73,13 +73,13 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return; return;
} }
static double dotVector3( double* vectorX, double* vectorY ){ static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2]; return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
} }
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();
...@@ -89,15 +89,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -89,15 +89,7 @@ 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 };
...@@ -108,7 +100,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -108,7 +100,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
// and various intermediate terms // and various intermediate terms
double deltaR[LastDeltaIndex][6]; double deltaR[LastDeltaIndex][6];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[AB][ii] = positions[particle1][ii] - positions[particle2][ii]; deltaR[AB][ii] = positions[particle1][ii] - positions[particle2][ii];
deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii]; deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii];
deltaR[DB][ii] = positions[particle4][ii] - positions[particle2][ii]; deltaR[DB][ii] = positions[particle4][ii] - positions[particle2][ii];
...@@ -116,38 +108,31 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -116,38 +108,31 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii]; deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii];
} }
double rDB2 = dotVector3( deltaR[DB], deltaR[DB] ); double rDB2 = dotVector3(deltaR[DB], deltaR[DB]);
double rAD2 = dotVector3( deltaR[AD], deltaR[AD] ); double rAD2 = dotVector3(deltaR[AD], deltaR[AD]);
double rCD2 = dotVector3( deltaR[CD], deltaR[CD] ); double rCD2 = dotVector3(deltaR[CD], deltaR[CD]);
double tempVector[3]; double tempVector[3];
crossProductVector3( deltaR[CB], deltaR[DB], tempVector ); crossProductVector3(deltaR[CB], deltaR[DB], tempVector);
double eE = dotVector3( deltaR[AB], tempVector ); double eE = dotVector3(deltaR[AB], tempVector );
double dot = dotVector3( deltaR[AD], deltaR[CD] ); double dot = dotVector3(deltaR[AD], deltaR[CD]);
double cc = rAD2*rCD2 - dot*dot; double cc = rAD2*rCD2 - dot*dot;
if( rDB2 <= 0.0 || cc == 0.0 ){ if (rDB2 <= 0.0 || cc == 0.0) {
return; return;
} }
double bkk2 = rDB2 - eE*eE/cc; double bkk2 = rDB2 - eE*eE/cc;
double cosine = sqrt(bkk2/rDB2); double cosine = sqrt(bkk2/rDB2);
double angle; double angle;
if( cosine >= 1.0 ){ if (cosine >= 1.0) {
angle = 0.0; angle = 0.0;
} }
else if( cosine <= -1.0 ){ else if (cosine <= -1.0) {
angle = PI_M; angle = PI_M;
} }
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
...@@ -161,23 +146,23 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -161,23 +146,23 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
double dEdCos; double dEdCos;
dEdCos = dEdDt/sqrt(cc*bkk2); dEdCos = dEdDt/sqrt(cc*bkk2);
if( eE > 0.0 ){ if (eE > 0.0) {
dEdCos *= -1.0; dEdCos *= -1.0;
} }
double term = eE/cc; double term = eE/cc;
double dccd[LastAtomIndex][3]; double dccd[LastAtomIndex][3];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
dccd[A][ii] = (deltaR[AD][ii]*rCD2 - deltaR[CD][ii]*dot)*term; dccd[A][ii] = (deltaR[AD][ii]*rCD2 - deltaR[CD][ii]*dot)*term;
dccd[C][ii] = (deltaR[CD][ii]*rAD2 - deltaR[AD][ii]*dot)*term; dccd[C][ii] = (deltaR[CD][ii]*rAD2 - deltaR[AD][ii]*dot)*term;
dccd[D][ii] = -1.0*(dccd[A][ii] + dccd[C][ii]); dccd[D][ii] = -1.0*(dccd[A][ii] + dccd[C][ii]);
} }
double deed[LastAtomIndex][3]; double deed[LastAtomIndex][3];
crossProductVector3( deltaR[DB], deltaR[CB], deed[A] ); crossProductVector3(deltaR[DB], deltaR[CB], deed[A]);
crossProductVector3( deltaR[AB], deltaR[DB], deed[C] ); crossProductVector3(deltaR[AB], deltaR[DB], deed[C]);
crossProductVector3( deltaR[CB], deltaR[AB], deed[D] ); crossProductVector3(deltaR[CB], deltaR[AB], deed[D]);
term = eE/rDB2; term = eE/rDB2;
deed[D][0] += deltaR[DB][0]*term; deed[D][0] += deltaR[DB][0]*term;
...@@ -189,24 +174,24 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -189,24 +174,24 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
// forces // forces
// calculate forces for atoms a, c, d // calculate forces for atoms a, c, d
// the force for b is then -( a+ c + d) // the force for b is then -(a+ c + d)
double subForce[LastAtomIndex][3]; double subForce[LastAtomIndex][3];
for( int jj = 0; jj < LastAtomIndex; jj++ ){ for (int jj = 0; jj < LastAtomIndex; jj++) {
// A, C, D // A, C, D
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
subForce[jj][ii] = dEdCos*( dccd[jj][ii] + deed[jj][ii] ); subForce[jj][ii] = dEdCos*(dccd[jj][ii] + deed[jj][ii]);
} }
if( jj == 0 )jj++; // skip B if (jj == 0)jj++; // skip B
// now compute B // now compute B
if( jj == 3 ){ if (jj == 3) {
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
subForce[1][ii] = -1.0*(subForce[0][ii] + subForce[2][ii] + subForce[3][ii]); subForce[1][ii] = -1.0*(subForce[0][ii] + subForce[2][ii] + subForce[3][ii]);
} }
} }
...@@ -243,70 +228,48 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -243,70 +228,48 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
return; return;
} }
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
State state = context.getState(State::Positions); State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions(); std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() ); expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){ for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0; expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
} }
// calculates forces/energy // calculates forces/energy
*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);
} }
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;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
...@@ -314,29 +277,29 @@ void testOneOutOfPlaneBend( FILE* log ) { ...@@ -314,29 +277,29 @@ void testOneOutOfPlaneBend( FILE* log ) {
AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce(); AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce();
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic( 0.5600000E-04 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic( 0.5600000E-04);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic( -0.7000000E-06 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic( -0.7000000E-06);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07);
double kOutOfPlaneBend = 0.328682196E-01; double kOutOfPlaneBend = 0.328682196E-01;
amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend ); amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend);
system.addForce(amoebaOutOfPlaneBendForce); system.addForce(amoebaOutOfPlaneBendForce);
ASSERT(!amoebaOutOfPlaneBendForce->usesPeriodicBoundaryConditions()); ASSERT(!amoebaOutOfPlaneBendForce->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions()); ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, Platform::getPlatformByName( "Reference")); Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 0.262660000E+02, 0.254130000E+02, 0.284200000E+01 ); positions[0] = Vec3(0.262660000E+02, 0.254130000E+02, 0.284200000E+01);
positions[1] = Vec3( 0.269130000E+02, 0.266390000E+02, 0.353100000E+01 ); positions[1] = Vec3(0.269130000E+02, 0.266390000E+02, 0.353100000E+01);
positions[2] = Vec3( 0.278860000E+02, 0.264630000E+02, 0.426300000E+01 ); positions[2] = Vec3(0.278860000E+02, 0.264630000E+02, 0.426300000E+01);
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,21 +307,21 @@ void testOneOutOfPlaneBend( FILE* log ) { ...@@ -344,21 +307,21 @@ 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;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
...@@ -366,10 +329,10 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -366,10 +329,10 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce(); AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce();
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic( 0.5600000E-04 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic( 0.5600000E-04);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic( -0.7000000E-06 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic( -0.7000000E-06);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07 ); amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07);
/* /*
285 441 442 443 444 0.328682196E-01 285 441 442 443 444 0.328682196E-01
286 441 442 444 443 0.164493407E-01 286 441 442 444 443 0.164493407E-01
...@@ -389,166 +352,125 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -389,166 +352,125 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
*/ */
std::map<int,Vec3> coordinates; std::map<int,Vec3> coordinates;
coordinates[440] = Vec3( 0.893800000E+01, 0.439800000E+01, 0.343100000E+01 ); coordinates[440] = Vec3( 0.893800000E+01, 0.439800000E+01, 0.343100000E+01);
coordinates[441] = Vec3( 0.779100000E+01, 0.614600000E+01, 0.390100000E+01 ); coordinates[441] = Vec3( 0.779100000E+01, 0.614600000E+01, 0.390100000E+01);
coordinates[442] = Vec3( 0.915400000E+01, 0.683900000E+01, 0.389400000E+01 ); coordinates[442] = Vec3( 0.915400000E+01, 0.683900000E+01, 0.389400000E+01);
coordinates[443] = Vec3( 0.101770000E+02, 0.619000000E+01, 0.379900000E+01 ); coordinates[443] = Vec3( 0.101770000E+02, 0.619000000E+01, 0.379900000E+01);
coordinates[444] = Vec3( 0.921000000E+01, 0.813800000E+01, 0.398600000E+01 ); coordinates[444] = Vec3( 0.921000000E+01, 0.813800000E+01, 0.398600000E+01);
coordinates[445] = Vec3( 0.708500000E+01, 0.672900000E+01, 0.332700000E+01 ); coordinates[445] = Vec3( 0.708500000E+01, 0.672900000E+01, 0.332700000E+01);
coordinates[446] = Vec3( 0.744300000E+01, 0.605200000E+01, 0.491900000E+01 ); coordinates[446] = Vec3( 0.744300000E+01, 0.605200000E+01, 0.491900000E+01);
coordinates[447] = Vec3( 0.100820000E+02, 0.859300000E+01, 0.398200000E+01 ); coordinates[447] = Vec3( 0.100820000E+02, 0.859300000E+01, 0.398200000E+01);
coordinates[448] = Vec3( 0.838000000E+01, 0.866100000E+01, 0.406000000E+01 ); coordinates[448] = Vec3( 0.838000000E+01, 0.866100000E+01, 0.406000000E+01);
double kOutOfPlaneBend = 0.328682196E-01; double kOutOfPlaneBend = 0.328682196E-01;
std::vector<int> particleIndices; std::vector<int> particleIndices;
if( setId == 1 ){ if (setId == 1) {
particleIndices.push_back( 441 ); particleIndices.push_back(441);
particleIndices.push_back( 442 ); particleIndices.push_back(442);
particleIndices.push_back( 443 ); particleIndices.push_back(443);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
kOutOfPlaneBend = 0.328682196E-01; kOutOfPlaneBend = 0.328682196E-01;
} }
else if( setId == 2 ){ else if (setId == 2) {
particleIndices.push_back( 441 ); particleIndices.push_back(441);
particleIndices.push_back( 442 ); particleIndices.push_back(442);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
particleIndices.push_back( 443 ); particleIndices.push_back(443);
kOutOfPlaneBend = 0.164493407E-01; kOutOfPlaneBend = 0.164493407E-01;
} }
else if( setId == 3 ){ else if (setId == 3) {
particleIndices.push_back( 443 ); particleIndices.push_back(443);
particleIndices.push_back( 442 ); particleIndices.push_back(442);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
particleIndices.push_back( 441 ); particleIndices.push_back(441);
kOutOfPlaneBend = 0.636650407E-02; kOutOfPlaneBend = 0.636650407E-02;
} }
else if( setId == 4 ){ else if (setId == 4) {
particleIndices.push_back( 442 ); particleIndices.push_back(442);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
particleIndices.push_back( 447 ); particleIndices.push_back(447);
particleIndices.push_back( 448 ); particleIndices.push_back(448);
kOutOfPlaneBend = 0.392956472E-02; kOutOfPlaneBend = 0.392956472E-02;
} }
else if( setId == 5 ){ else if (setId == 5) {
particleIndices.push_back( 442 ); particleIndices.push_back(442);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
particleIndices.push_back( 448 ); particleIndices.push_back(448);
particleIndices.push_back( 447 ); particleIndices.push_back(447);
kOutOfPlaneBend = 0.392956472E-02; kOutOfPlaneBend = 0.392956472E-02;
} }
else if( setId == 6 ){ else if (setId == 6) {
particleIndices.push_back( 447 ); particleIndices.push_back(447);
particleIndices.push_back( 444 ); particleIndices.push_back(444);
particleIndices.push_back( 448 ); particleIndices.push_back(448);
particleIndices.push_back( 442 ); particleIndices.push_back(442);
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);
system.addForce(amoebaOutOfPlaneBendForce); system.addForce(amoebaOutOfPlaneBendForce);
Context context(system, integrator, Platform::getPlatformByName( "Reference")); Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
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());
} }
positions[ii] = coordinates[particleIndices[ii]]; positions[ii] = coordinates[particleIndices[ii]];
} }
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;
static double totalEnergy; static double totalEnergy;
if( iter == 0 ){ if (iter == 0) {
totalForces[441] = Vec3( 0.0, 0.0, 0.0 ); totalForces[441] = Vec3( 0.0, 0.0, 0.0);
totalForces[442] = Vec3( 0.0, 0.0, 0.0 ); totalForces[442] = Vec3( 0.0, 0.0, 0.0);
totalForces[443] = Vec3( 0.0, 0.0, 0.0 ); totalForces[443] = Vec3( 0.0, 0.0, 0.0);
totalForces[444] = Vec3( 0.0, 0.0, 0.0 ); totalForces[444] = Vec3( 0.0, 0.0, 0.0);
totalForces[445] = Vec3( 0.0, 0.0, 0.0 ); totalForces[445] = Vec3( 0.0, 0.0, 0.0);
totalForces[446] = Vec3( 0.0, 0.0, 0.0 ); totalForces[446] = Vec3( 0.0, 0.0, 0.0);
totalForces[447] = Vec3( 0.0, 0.0, 0.0 ); totalForces[447] = Vec3( 0.0, 0.0, 0.0);
totalForces[448] = Vec3( 0.0, 0.0, 0.0 ); totalForces[448] = Vec3( 0.0, 0.0, 0.0);
totalForces[449] = Vec3( 0.0, 0.0, 0.0 ); totalForces[449] = Vec3( 0.0, 0.0, 0.0);
totalEnergy = 0.0; totalEnergy = 0.0;
} }
iter++; iter++;
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++) {
for( unsigned int jj = 0; jj < 3; jj++ ){ for (unsigned int jj = 0; jj < 3; jj++) {
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[]) {
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; //for (int ii = 1; ii <= 6; ii++) {
//FILE* log = fopen( "AmoebaOutOfPlaneBendForce.log", "w" );; // testOneOutOfPlaneBend2(ii);
testOneOutOfPlaneBend( log );
//testOneOutOfPlaneBend2( log, atoi( argv[1] ) );
//for( int ii = 1; ii <= 6; ii++ ){
// testOneOutOfPlaneBend2( log, 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 @ 162d69a2
...@@ -63,7 +63,7 @@ const double TOL = 1e-5; ...@@ -63,7 +63,7 @@ const double TOL = 1e-5;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){ static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1]; vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2]; vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
...@@ -72,25 +72,18 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto ...@@ -72,25 +72,18 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return; return;
} }
static double dotVector3( double* vectorX, double* vectorY ){ static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2]; return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
} }
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];
...@@ -98,16 +91,16 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit ...@@ -98,16 +91,16 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
enum { A, B, C, D, E, F, LastAtomIndex }; enum { A, B, C, D, E, F, LastAtomIndex };
double d[LastAtomIndex][3]; double d[LastAtomIndex][3];
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii]; deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii];
deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii]; deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii];
deltaR[EC][ii] = positions[particle5][ii] - positions[particle3][ii]; deltaR[EC][ii] = positions[particle5][ii] - positions[particle3][ii];
deltaR[FC][ii] = positions[particle6][ii] - positions[particle3][ii]; deltaR[FC][ii] = positions[particle6][ii] - positions[particle3][ii];
} }
crossProductVector3( deltaR[AD], deltaR[BD], deltaR[P] ); crossProductVector3(deltaR[AD], deltaR[BD], deltaR[P]);
crossProductVector3( deltaR[EC], deltaR[FC], deltaR[Q] ); crossProductVector3(deltaR[EC], deltaR[FC], deltaR[Q]);
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[CP][ii] = -deltaR[P][ii]; deltaR[CP][ii] = -deltaR[P][ii];
deltaR[DC][ii] = positions[particle4][ii] - positions[particle3][ii]; deltaR[DC][ii] = positions[particle4][ii] - positions[particle3][ii];
deltaR[QD][ii] = deltaR[Q][ii]; deltaR[QD][ii] = deltaR[Q][ii];
...@@ -115,25 +108,25 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit ...@@ -115,25 +108,25 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
deltaR[P][ii] += positions[particle3][ii]; deltaR[P][ii] += positions[particle3][ii];
deltaR[Q][ii] += positions[particle4][ii]; deltaR[Q][ii] += positions[particle4][ii];
} }
crossProductVector3( deltaR[CP], deltaR[DC], deltaR[T] ); crossProductVector3(deltaR[CP], deltaR[DC], deltaR[T] );
crossProductVector3( deltaR[DC], deltaR[QD], deltaR[U] ); crossProductVector3(deltaR[DC], deltaR[QD], deltaR[U] );
crossProductVector3( deltaR[T], deltaR[U], deltaR[TU] ); crossProductVector3(deltaR[T], deltaR[U], deltaR[TU]);
double rT2 = dotVector3( deltaR[T], deltaR[T] ); double rT2 = dotVector3(deltaR[T], deltaR[T]);
double rU2 = dotVector3( deltaR[U], deltaR[U] ); double rU2 = dotVector3(deltaR[U], deltaR[U]);
double rTrU = sqrt( rT2*rU2 ); double rTrU = sqrt(rT2*rU2);
if( rTrU <= 0.0 ){ if (rTrU <= 0.0) {
return; return;
} }
double rDC = dotVector3( deltaR[DC], deltaR[DC] ); double rDC = dotVector3(deltaR[DC], deltaR[DC]);
rDC = sqrt( rDC ); rDC = sqrt(rDC);
double cosine = dotVector3( deltaR[T], deltaR[U] ); double cosine = dotVector3(deltaR[T], deltaR[U]);
cosine /= rTrU; cosine /= rTrU;
double sine = dotVector3( deltaR[DC], deltaR[TU] ); double sine = dotVector3(deltaR[DC], deltaR[TU]);
sine /= ( rDC*rTrU ); sine /= (rDC*rTrU);
double cosine2 = cosine*cosine - sine*sine; double cosine2 = cosine*cosine - sine*sine;
double sine2 = 2.0*cosine*sine; double sine2 = 2.0*cosine*sine;
...@@ -143,37 +136,37 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit ...@@ -143,37 +136,37 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
double dedphi = kTorsion*dphi2; double dedphi = kTorsion*dphi2;
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[DP][ii] = positions[particle4][ii] - deltaR[P][ii]; deltaR[DP][ii] = positions[particle4][ii] - deltaR[P][ii];
deltaR[QC][ii] = deltaR[Q][ii] - positions[particle3][ii]; deltaR[QC][ii] = deltaR[Q][ii] - positions[particle3][ii];
} }
double factorT = dedphi/( rDC*rT2 ); double factorT = dedphi/(rDC*rT2);
double factorU = -dedphi/( rDC*rU2 ); double factorU = -dedphi/(rDC*rU2);
crossProductVector3( deltaR[T], deltaR[DC], deltaR[dT] ); crossProductVector3(deltaR[T], deltaR[DC], deltaR[dT] );
crossProductVector3( deltaR[U], deltaR[DC], deltaR[dU] ); crossProductVector3(deltaR[U], deltaR[DC], deltaR[dU] );
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[dT][ii] *= factorT; deltaR[dT][ii] *= factorT;
deltaR[dU][ii] *= factorU; deltaR[dU][ii] *= factorU;
} }
crossProductVector3( deltaR[dT], deltaR[DC], deltaR[dP] ); crossProductVector3(deltaR[dT], deltaR[DC], deltaR[dP] );
crossProductVector3( deltaR[dU], deltaR[DC], deltaR[dQ] ); crossProductVector3(deltaR[dU], deltaR[DC], deltaR[dQ] );
crossProductVector3( deltaR[DP], deltaR[dT], deltaR[dC1] ); crossProductVector3(deltaR[DP], deltaR[dT], deltaR[dC1] );
crossProductVector3( deltaR[dU], deltaR[QD], deltaR[dC2] ); crossProductVector3(deltaR[dU], deltaR[QD], deltaR[dC2] );
crossProductVector3( deltaR[dT], deltaR[CP], deltaR[dD1] ); crossProductVector3(deltaR[dT], deltaR[CP], deltaR[dD1] );
crossProductVector3( deltaR[QC], deltaR[dU], deltaR[dD2] ); crossProductVector3(deltaR[QC], deltaR[dU], deltaR[dD2] );
crossProductVector3( deltaR[BD], deltaR[dP], d[A] ); crossProductVector3(deltaR[BD], deltaR[dP], d[A] );
crossProductVector3( deltaR[dP], deltaR[AD], d[B] ); crossProductVector3(deltaR[dP], deltaR[AD], d[B] );
crossProductVector3( deltaR[FC], deltaR[dQ], d[E] ); crossProductVector3(deltaR[FC], deltaR[dQ], d[E] );
crossProductVector3( deltaR[dQ], deltaR[EC], d[F] ); crossProductVector3(deltaR[dQ], deltaR[EC], d[F] );
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
d[C][ii] = deltaR[dC1][ii] + deltaR[dC2][ii] + deltaR[dP][ii] - d[E][ii] - d[F][ii]; d[C][ii] = deltaR[dC1][ii] + deltaR[dC2][ii] + deltaR[dP][ii] - d[E][ii] - d[F][ii];
d[D][ii] = deltaR[dD1][ii] + deltaR[dD2][ii] + deltaR[dQ][ii] - d[A][ii] - d[B][ii]; d[D][ii] = deltaR[dD1][ii] + deltaR[dD2][ii] + deltaR[dQ][ii] - d[A][ii] - d[B][ii];
} }
...@@ -211,69 +204,48 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit ...@@ -211,69 +204,48 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
return; return;
} }
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
State state = context.getState(State::Positions); State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions(); std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() ); expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){ for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0; expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
} }
// calculates forces/energy // calculates forces/energy
*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);
} }
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;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
...@@ -282,25 +254,25 @@ void testOnePiTorsion( FILE* log ) { ...@@ -282,25 +254,25 @@ void testOnePiTorsion( FILE* log ) {
AmoebaPiTorsionForce* amoebaPiTorsionForce = new AmoebaPiTorsionForce(); AmoebaPiTorsionForce* amoebaPiTorsionForce = new AmoebaPiTorsionForce();
double kTorsion = 6.85; double kTorsion = 6.85;
amoebaPiTorsionForce->addPiTorsion(0, 1, 2, 3, 4, 5, kTorsion ); amoebaPiTorsionForce->addPiTorsion(0, 1, 2, 3, 4, 5, kTorsion);
system.addForce(amoebaPiTorsionForce); system.addForce(amoebaPiTorsionForce);
ASSERT(!amoebaPiTorsionForce->usesPeriodicBoundaryConditions()); ASSERT(!amoebaPiTorsionForce->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions()); ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, Platform::getPlatformByName( "Reference")); Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 0.262660000E+02, 0.254130000E+02, 0.284200000E+01 ); positions[0] = Vec3(0.262660000E+02, 0.254130000E+02, 0.284200000E+01);
positions[1] = Vec3( 0.278860000E+02, 0.264630000E+02, 0.426300000E+01 ); positions[1] = Vec3(0.278860000E+02, 0.264630000E+02, 0.426300000E+01);
positions[2] = Vec3( 0.269130000E+02, 0.266390000E+02, 0.353100000E+01 ); positions[2] = Vec3(0.269130000E+02, 0.266390000E+02, 0.353100000E+01);
positions[3] = Vec3( 0.245568230E+02, 0.250215290E+02, 0.796852800E+01 ); positions[3] = Vec3(0.245568230E+02, 0.250215290E+02, 0.796852800E+01);
positions[4] = Vec3( 0.261000000E+02, 0.292530000E+02, 0.520200000E+01 ); positions[4] = Vec3(0.261000000E+02, 0.292530000E+02, 0.520200000E+01);
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,31 +280,22 @@ void testOnePiTorsion( FILE* log ) { ...@@ -308,31 +280,22 @@ 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[]) {
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 @ 162d69a2
...@@ -65,7 +65,7 @@ const double TOL = 1e-4; ...@@ -65,7 +65,7 @@ const double TOL = 1e-4;
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){ static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1]; vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2]; vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
...@@ -74,26 +74,19 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto ...@@ -74,26 +74,19 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return; return;
} }
static double dotVector3( double* vectorX, double* vectorY ){ static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2]; return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
} }
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 };
...@@ -106,31 +99,31 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -106,31 +99,31 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
double deltaR[LastDeltaIndex][3]; double deltaR[LastDeltaIndex][3];
double rAB2 = 0.0; double rAB2 = 0.0;
double rCB2 = 0.0; double rCB2 = 0.0;
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[AB][ii] = positions[particle1][ii] - positions[particle2][ii]; deltaR[AB][ii] = positions[particle1][ii] - positions[particle2][ii];
rAB2 += deltaR[AB][ii]*deltaR[AB][ii]; rAB2 += deltaR[AB][ii]*deltaR[AB][ii];
deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii]; deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii];
rCB2 += deltaR[CB][ii]*deltaR[CB][ii]; rCB2 += deltaR[CB][ii]*deltaR[CB][ii];
} }
double rAB = sqrt( rAB2 ); double rAB = sqrt(rAB2);
double rCB = sqrt( rCB2 ); double rCB = sqrt(rCB2);
crossProductVector3( deltaR[CB], deltaR[AB], deltaR[CBxAB] ); crossProductVector3(deltaR[CB], deltaR[AB], deltaR[CBxAB]);
double rP = dotVector3( deltaR[CBxAB], deltaR[CBxAB] ); double rP = dotVector3(deltaR[CBxAB], deltaR[CBxAB]);
rP = sqrt( rP ); rP = sqrt(rP);
if( rP <= 0.0 ){ if (rP <= 0.0) {
return; return;
} }
double dot = dotVector3( deltaR[CB], deltaR[AB] ); double dot = dotVector3(deltaR[CB], deltaR[AB]);
double cosine = dot/(rAB*rCB); double cosine = dot/(rAB*rCB);
double angle; double angle;
if( cosine >= 1.0 ){ if (cosine >= 1.0) {
angle = 0.0; angle = 0.0;
} }
else if( cosine <= -1.0 ){ else if (cosine <= -1.0) {
angle = PI_M; angle = PI_M;
} }
else { else {
...@@ -142,9 +135,9 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -142,9 +135,9 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
// P = CBxAB // P = CBxAB
crossProductVector3( deltaR[AB], deltaR[CBxAB], deltaR[ABxP] ); crossProductVector3(deltaR[AB], deltaR[CBxAB], deltaR[ABxP]);
crossProductVector3( deltaR[CB], deltaR[CBxAB], deltaR[CBxP] ); crossProductVector3(deltaR[CB], deltaR[CBxAB], deltaR[CBxP]);
for( int ii = 0; ii < 3; ii++ ){ for (int ii = 0; ii < 3; ii++) {
deltaR[ABxP][ii] *= termA; deltaR[ABxP][ii] *= termA;
deltaR[CBxP][ii] *= termC; deltaR[CBxP][ii] *= termC;
} }
...@@ -162,14 +155,14 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -162,14 +155,14 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
// forces // forces
// calculate forces for atoms a, b, c // calculate forces for atoms a, b, c
// the force for b is then -( a + c) // the force for b is then -(a + c)
double subForce[LastAtomIndex][3]; double subForce[LastAtomIndex][3];
double dt = angle - angleStretchBend; double dt = angle - angleStretchBend;
for( int jj = 0; jj < 3; jj++ ){ for (int jj = 0; jj < 3; jj++) {
subForce[A][jj] = kStretchBend*dt*termA*deltaR[AB][jj] + drkk*deltaR[ABxP][jj]; subForce[A][jj] = kStretchBend*dt*termA*deltaR[AB][jj] + drkk*deltaR[ABxP][jj];
subForce[C][jj] = k2StretchBend*dt*termC*deltaR[CB][jj] + drkk*deltaR[CBxP][jj]; subForce[C][jj] = k2StretchBend*dt*termC*deltaR[CB][jj] + drkk*deltaR[CBxP][jj];
subForce[B][jj] = -( subForce[A][jj] + subForce[C][jj] ); subForce[B][jj] = -(subForce[A][jj] + subForce[C][jj]);
} }
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -189,79 +182,49 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -189,79 +182,49 @@ 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
State state = context.getState(State::Positions); State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions(); std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() ); expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){ for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0; expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
} }
// calculates forces/energy // calculates forces/energy
*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 for (unsigned int ii = 0; ii < forces.size(); ii++) {
if( log ){ ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
(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++ ){
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;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
...@@ -275,21 +238,21 @@ void testOneStretchBend( FILE* log ) { ...@@ -275,21 +238,21 @@ void testOneStretchBend( FILE* log ) {
//double kStretchBend = 0.750491578E-01; //double kStretchBend = 0.750491578E-01;
double kStretchBend = 1.0; double kStretchBend = 1.0;
amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend ); amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend);
system.addForce(amoebaStretchBendForce); system.addForce(amoebaStretchBendForce);
ASSERT(!amoebaStretchBendForce->usesPeriodicBoundaryConditions()); ASSERT(!amoebaStretchBendForce->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions()); ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, Platform::getPlatformByName( "Reference")); Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 0.262660000E+02, 0.254130000E+02, 0.284200000E+01 ); positions[0] = Vec3(0.262660000E+02, 0.254130000E+02, 0.284200000E+01);
positions[1] = Vec3( 0.273400000E+02, 0.244300000E+02, 0.261400000E+01 ); positions[1] = Vec3(0.273400000E+02, 0.244300000E+02, 0.261400000E+01);
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,31 +260,22 @@ void testOneStretchBend( FILE* log ) { ...@@ -297,31 +260,22 @@ 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[]) {
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/TestReferenceAmoebaTorsionTorsionForce.cpp
View file @ 162d69a2
...@@ -2559,20 +2559,20 @@ static double grid[4][625][6] = { ...@@ -2559,20 +2559,20 @@ static double grid[4][625][6] = {
int elementCount = (includeDerivs ? 6 : 3); int elementCount = (includeDerivs ? 6 : 3);
std::vector<TorsionTorsionGrid> grids(4); std::vector<TorsionTorsionGrid> grids(4);
for( int ii = 0; ii < 4; ii++ ){ for (int ii = 0; ii < 4; ii++) {
grids[ii].resize( 25 ); grids[ii].resize(25);
for( int jj = 0; jj < 25; jj++ ){ for (int jj = 0; jj < 25; jj++) {
grids[ii][jj].resize(25); grids[ii][jj].resize(25);
for( int kk = 0; kk < 25; kk++ ){ for (int kk = 0; kk < 25; kk++) {
grids[ii][jj][kk].resize(elementCount); grids[ii][jj][kk].resize(elementCount);
} }
} }
int index = 0; int index = 0;
for( int jj = 0; jj < 25; jj++ ){ for (int jj = 0; jj < 25; jj++) {
for( int kk = 0; kk < 25; kk++ ){ for (int kk = 0; kk < 25; kk++) {
int jjIndex = static_cast<int>(((grid[ii][index][0] + 180.0)/15.0)+1.0e-05); int jjIndex = static_cast<int>(((grid[ii][index][0] + 180.0)/15.0)+1.0e-05);
int kkIndex = static_cast<int>(((grid[ii][index][1] + 180.0)/15.0)+1.0e-05); int kkIndex = static_cast<int>(((grid[ii][index][1] + 180.0)/15.0)+1.0e-05);
for( int ll = 0; ll < elementCount; ll++ ){ for (int ll = 0; ll < elementCount; ll++) {
grids[ii][kk][jj][ll] = grid[ii][index][ll]; grids[ii][kk][jj][ll] = grid[ii][index][ll];
} }
index++; index++;
...@@ -2586,7 +2586,7 @@ void testTorsionTorsion(int systemId, bool includeDerivs) { ...@@ -2586,7 +2586,7 @@ void testTorsionTorsion(int systemId, bool includeDerivs) {
System system; System system;
int numberOfParticles = 6; int numberOfParticles = 6;
for( int ii = 0; ii < numberOfParticles; ii++ ){ for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0); system.addParticle(1.0);
} }
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
...@@ -2598,48 +2598,48 @@ void testTorsionTorsion(int systemId, bool includeDerivs) { ...@@ -2598,48 +2598,48 @@ void testTorsionTorsion(int systemId, bool includeDerivs) {
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy; double expectedEnergy;
if( systemId == 0 ){ if (systemId == 0) {
// villin: 2 19 20 21 38 25 grid=2 // villin: 2 19 20 21 38 25 grid=2
chiralCheckAtomIndex = 5; chiralCheckAtomIndex = 5;
gridIndex = 2; gridIndex = 2;
positions[0] = Vec3( -0.422792800E+01, -0.110605910E+02, -0.508156700E+01 ); positions[0] = Vec3(-0.422792800E+01, -0.110605910E+02, -0.508156700E+01);
positions[1] = Vec3( -0.447153100E+01, -0.978627900E+01, -0.466405800E+01 ); positions[1] = Vec3(-0.447153100E+01, -0.978627900E+01, -0.466405800E+01);
positions[2] = Vec3( -0.531878400E+01, -0.940508600E+01, -0.352283100E+01 ); positions[2] = Vec3(-0.531878400E+01, -0.940508600E+01, -0.352283100E+01);
positions[3] = Vec3( -0.679606000E+01, -0.974353100E+01, -0.382975700E+01 ); positions[3] = Vec3(-0.679606000E+01, -0.974353100E+01, -0.382975700E+01);
positions[4] = Vec3( -0.760612300E+01, -0.992590200E+01, -0.275088400E+01 ); positions[4] = Vec3(-0.760612300E+01, -0.992590200E+01, -0.275088400E+01);
positions[5] = Vec3( -0.516893900E+01, -0.788347000E+01, -0.316943000E+01 ); positions[5] = Vec3(-0.516893900E+01, -0.788347000E+01, -0.316943000E+01);
expectedForces[0] = Vec3( 0.906091624E+00, -0.529814945E-01, 0.690384140E+00 ); expectedForces[0] = Vec3( 0.906091624E+00, -0.529814945E-01, 0.690384140E+00);
expectedForces[1] = Vec3( -0.124550232E+01, -0.999341692E+00, -0.590867130E+00 ); expectedForces[1] = Vec3(-0.124550232E+01, -0.999341692E+00, -0.590867130E+00);
expectedForces[2] = Vec3( 0.534419689E+00, 0.612404926E-01, 0.547380310E-01 ); expectedForces[2] = Vec3( 0.534419689E+00, 0.612404926E-01, 0.547380310E-01);
expectedForces[3] = Vec3( -5.732010432E-01, 2.645718463E+00, -1.585204274E-01 ); expectedForces[3] = Vec3(-5.732010432E-01, 2.645718463E+00, -1.585204274E-01);
expectedForces[4] = Vec3( 3.781920539E-01, -1.654635768E+00, 4.265386268E-03 ); expectedForces[4] = Vec3( 3.781920539E-01, -1.654635768E+00, 4.265386268E-03);
expectedForces[5] = Vec3( 0.0, 0.0, 0.0 ); expectedForces[5] = Vec3( 0.0, 0.0, 0.0);
expectedEnergy = -2.699654759E+00; expectedEnergy = -2.699654759E+00;
} }
else if( systemId == 1 ){ else if (systemId == 1) {
// villin: 158 176 177 178 183 -1 0 // villin: 158 176 177 178 183 -1 0
chiralCheckAtomIndex = -1; chiralCheckAtomIndex = -1;
gridIndex = 0; gridIndex = 0;
positions[0] = Vec3( -0.105946640E+02, -0.917797000E+00, 0.105486310E+02 ); positions[0] = Vec3(-0.105946640E+02, -0.917797000E+00, 0.105486310E+02);
positions[1] = Vec3( -0.115059090E+02, -0.141876700E+01, 0.966933200E+01 ); positions[1] = Vec3(-0.115059090E+02, -0.141876700E+01, 0.966933200E+01);
positions[2] = Vec3( -0.128314660E+02, -0.876338000E+00, 0.942959800E+01 ); positions[2] = Vec3(-0.128314660E+02, -0.876338000E+00, 0.942959800E+01);
positions[3] = Vec3( -0.130879850E+02, -0.760280000E-01, 0.814732200E+01 ); positions[3] = Vec3(-0.130879850E+02, -0.760280000E-01, 0.814732200E+01);
positions[4] = Vec3( -0.120888080E+02, 0.112050000E-01, 0.722704500E+01 ); positions[4] = Vec3(-0.120888080E+02, 0.112050000E-01, 0.722704500E+01);
positions[5] = Vec3( 0.0, 0.0, 0.0 ); positions[5] = Vec3( 0.0, 0.0, 0.0);
expectedForces[0] = Vec3( 4.165851130E-01, 6.608242922E-01, -8.082168261E-01 ); expectedForces[0] = Vec3( 4.165851130E-01, 6.608242922E-01, -8.082168261E-01);
expectedForces[1] = Vec3( -6.024659721E-01, -8.878744406E-01, 1.322274444E+00 ); expectedForces[1] = Vec3(-6.024659721E-01, -8.878744406E-01, 1.322274444E+00);
expectedForces[2] = Vec3( 3.196925118E-02, -3.137497848E-01, -8.207984001E-01 ); expectedForces[2] = Vec3( 3.196925118E-02, -3.137497848E-01, -8.207984001E-01);
expectedForces[3] = Vec3( 3.842205941E-02, 2.602732089E-01, 1.547586195E-01 ); expectedForces[3] = Vec3( 3.842205941E-02, 2.602732089E-01, 1.547586195E-01);
expectedForces[4] = Vec3( 1.154895485E-01, 2.805267242E-01, 1.519821623E-01 ); expectedForces[4] = Vec3( 1.154895485E-01, 2.805267242E-01, 1.519821623E-01);
expectedForces[5] = Vec3( 0.0, 0.0, 0.0 ); expectedForces[5] = Vec3( 0.0, 0.0, 0.0);
expectedEnergy = -3.372536909E+00; expectedEnergy = -3.372536909E+00;
} }
...@@ -2656,21 +2656,21 @@ void testTorsionTorsion(int systemId, bool includeDerivs) { ...@@ -2656,21 +2656,21 @@ void testTorsionTorsion(int systemId, bool includeDerivs) {
std::vector<Vec3> forces = state.getForces(); std::vector<Vec3> forces = state.getForces();
const double conversion = -1.0; const double conversion = -1.0;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for (unsigned int ii = 0; ii < forces.size(); ii++) {
forces[ii][0] *= conversion; forces[ii][0] *= conversion;
forces[ii][1] *= conversion; forces[ii][1] *= conversion;
forces[ii][2] *= conversion; forces[ii][2] *= conversion;
} }
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);
} }
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance ); ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
} }
int main( int numberOfArguments, char* argv[] ) { int main(int numberOfArguments, char* argv[]) {
try { try {
std::cout << "TestReferenceAmoebaTorsionTorsionForce running test..." << std::endl; std::cout << "TestReferenceAmoebaTorsionTorsionForce running test..." << std::endl;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaVdwForce.cpp
View file @ 162d69a2
This source diff could not be displayed because it is too large. You can view the blob instead.
plugins/amoeba/platforms/reference/tests/TestReferenceWcaDispersionForce.cpp
View file @ 162d69a2
...@@ -55,21 +55,21 @@ extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories(); ...@@ -55,21 +55,21 @@ extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();
const double TOL = 1e-4; const double TOL = 1e-4;
void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy, 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, FILE* log) {
// beginning of WcaDispersion setup // beginning of WcaDispersion setup
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);
} }
ASSERT_EQUAL_TOL_MOD( expectedEnergy, energy, tolerance, testName ); ASSERT_EQUAL_TOL_MOD(expectedEnergy, energy, tolerance, testName);
} }
// test Wca dispersion // test Wca dispersion
void testWcaDispersionAmmonia( FILE* log ) { void testWcaDispersionAmmonia(FILE* log) {
std::string testName = "testWcaDispersionAmmonia"; std::string testName = "testWcaDispersionAmmonia";
int numberOfParticles = 8; int numberOfParticles = 8;
...@@ -79,41 +79,41 @@ void testWcaDispersionAmmonia( FILE* log ) { ...@@ -79,41 +79,41 @@ void testWcaDispersionAmmonia( FILE* log ) {
System system; System system;
AmoebaWcaDispersionForce* amoebaWcaDispersionForce = new AmoebaWcaDispersionForce();; AmoebaWcaDispersionForce* amoebaWcaDispersionForce = new AmoebaWcaDispersionForce();;
amoebaWcaDispersionForce->setEpso( 4.6024000e-01 ); amoebaWcaDispersionForce->setEpso( 4.6024000e-01);
amoebaWcaDispersionForce->setEpsh( 5.6484000e-02 ); amoebaWcaDispersionForce->setEpsh( 5.6484000e-02);
amoebaWcaDispersionForce->setRmino( 1.7025000e-01 ); amoebaWcaDispersionForce->setRmino( 1.7025000e-01);
amoebaWcaDispersionForce->setRminh( 1.3275000e-01 ); amoebaWcaDispersionForce->setRminh( 1.3275000e-01);
amoebaWcaDispersionForce->setDispoff( 2.6000000e-02 ); amoebaWcaDispersionForce->setDispoff(2.6000000e-02);
amoebaWcaDispersionForce->setAwater( 3.3428000e+01 ); amoebaWcaDispersionForce->setAwater( 3.3428000e+01);
amoebaWcaDispersionForce->setSlevy( 1.0000000e+00 ); amoebaWcaDispersionForce->setSlevy( 1.0000000e+00);
amoebaWcaDispersionForce->setShctd( 8.1000000e-01 ); amoebaWcaDispersionForce->setShctd( 8.1000000e-01);
// addParticle: radius, epsilon // addParticle: radius, epsilon
for( unsigned int ii = 0; ii < 2; ii++ ){ for (unsigned int ii = 0; ii < 2; ii++) {
system.addParticle( 1.4007000e+01 ); system.addParticle( 1.4007000e+01);
amoebaWcaDispersionForce->addParticle( 1.8550000e-01, 4.3932000e-01 ); amoebaWcaDispersionForce->addParticle( 1.8550000e-01, 4.3932000e-01);
system.addParticle( 1.0080000e+00 ); system.addParticle( 1.0080000e+00);
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02 ); amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02);
system.addParticle( 1.0080000e+00 ); system.addParticle( 1.0080000e+00);
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02 ); amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02);
system.addParticle( 1.0080000e+00 ); system.addParticle( 1.0080000e+00);
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02 ); amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02);
} }
std::vector<Vec3> positions(numberOfParticles); std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03 ); positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03);
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02 ); positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02);
positions[2] = Vec3( 2.0843610e-01, 8.0953200e-02, 3.7462200e-02 ); positions[2] = Vec3( 2.0843610e-01, 8.0953200e-02, 3.7462200e-02);
positions[3] = Vec3( 1.7280780e-01, 2.0730000e-04, -9.8741700e-02 ); positions[3] = Vec3( 1.7280780e-01, 2.0730000e-04, -9.8741700e-02);
positions[4] = Vec3( -1.6743680e-01, 1.5900000e-05, -6.6149000e-03 ); positions[4] = Vec3( -1.6743680e-01, 1.5900000e-05, -6.6149000e-03);
positions[5] = Vec3( -2.0428260e-01, 8.1071500e-02, 4.1343900e-02 ); positions[5] = Vec3( -2.0428260e-01, 8.1071500e-02, 4.1343900e-02);
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02 ); positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02);
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02 ); positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02);
system.addForce(amoebaWcaDispersionForce); system.addForce(amoebaWcaDispersionForce);
ASSERT(!amoebaWcaDispersionForce->usesPeriodicBoundaryConditions()); ASSERT(!amoebaWcaDispersionForce->usesPeriodicBoundaryConditions());
...@@ -122,7 +122,7 @@ void testWcaDispersionAmmonia( FILE* log ) { ...@@ -122,7 +122,7 @@ void testWcaDispersionAmmonia( FILE* log ) {
std::string platformName; std::string platformName;
platformName = "Reference"; platformName = "Reference";
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName( platformName ) ); Context context(system, integrator, Platform::getPlatformByName(platformName));
context.setPositions(positions); context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
...@@ -134,17 +134,17 @@ void testWcaDispersionAmmonia( FILE* log ) { ...@@ -134,17 +134,17 @@ void testWcaDispersionAmmonia( FILE* log ) {
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = -2.6981209e+01; double expectedEnergy = -2.6981209e+01;
expectedForces[0] = Vec3( 4.7839388e+00, -7.3510133e-04, -5.0382764e-01 ); expectedForces[0] = Vec3( 4.7839388e+00, -7.3510133e-04, -5.0382764e-01);
expectedForces[1] = Vec3( 1.4657758e+00, 1.2431003e+00, -6.7075886e-01 ); expectedForces[1] = Vec3( 1.4657758e+00, 1.2431003e+00, -6.7075886e-01);
expectedForces[2] = Vec3( 1.4563936e+00, -1.2399917e+00, -6.7443841e-01 ); expectedForces[2] = Vec3( 1.4563936e+00, -1.2399917e+00, -6.7443841e-01);
expectedForces[3] = Vec3( 2.1116744e+00, -2.7407512e-03, 1.3271245e+00 ); expectedForces[3] = Vec3( 2.1116744e+00, -2.7407512e-03, 1.3271245e+00);
expectedForces[4] = Vec3( -4.7528440e+00, -1.5148066e-03, 1.2653813e+00 ); expectedForces[4] = Vec3( -4.7528440e+00, -1.5148066e-03, 1.2653813e+00);
expectedForces[5] = Vec3( -1.1875619e+00, -1.2866678e+00, -3.9109060e-01 ); expectedForces[5] = Vec3( -1.1875619e+00, -1.2866678e+00, -3.9109060e-01);
expectedForces[6] = Vec3( -2.6885679e+00, -4.3038639e-04, 3.3763583e-02 ); expectedForces[6] = Vec3( -2.6885679e+00, -4.3038639e-04, 3.3763583e-02);
expectedForces[7] = Vec3( -1.1888087e+00, 1.2889802e+00, -3.8615387e-01 ); expectedForces[7] = Vec3( -1.1888087e+00, 1.2889802e+00, -3.8615387e-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, log);
// Try changing the particle parameters and make sure it's still correct. // Try changing the particle parameters and make sure it's still correct.
...@@ -172,7 +172,7 @@ void testWcaDispersionAmmonia( FILE* log ) { ...@@ -172,7 +172,7 @@ void testWcaDispersionAmmonia( FILE* log ) {
compareForcesEnergy(testName, state1.getPotentialEnergy(), state2.getPotentialEnergy(), state1.getForces(), state2.getForces(), tolerance, log); compareForcesEnergy(testName, state1.getPotentialEnergy(), state2.getPotentialEnergy(), state1.getForces(), state2.getForces(), tolerance, log);
} }
int main( int numberOfArguments, char* argv[] ) { int main(int numberOfArguments, char* argv[]) {
try { try {
std::cout << "TestCudaAmoebaWcaDispersionForce running test..." << std::endl; std::cout << "TestCudaAmoebaWcaDispersionForce running test..." << std::endl;
...@@ -182,7 +182,7 @@ int main( int numberOfArguments, char* argv[] ) { ...@@ -182,7 +182,7 @@ int main( int numberOfArguments, char* argv[] ) {
// test Wca dispersion force using two ammonia molecules // test Wca dispersion force using two ammonia molecules
testWcaDispersionAmmonia( log ); testWcaDispersionAmmonia(log);
} }
......
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