Skip to content

GitLab

Commit a568bb12 authored by peastman's avatar peastman
Browse files

Deleted more debugging code

parent 162d69a2
Show whitespace changes
Inline Side-by-side
Showing with 179 additions and 705 deletions
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaAngleForce.cpp
View file @ a568bb12
...@@ -77,7 +77,7 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto ...@@ -77,7 +77,7 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, double quadraticK, double cubicK, static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, 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 ){
...@@ -88,13 +88,6 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou ...@@ -88,13 +88,6 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou
angle = RADIAN*acos(cosine); angle = RADIAN*acos(cosine);
} }
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "getPrefactorsGivenAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealAngle );
(void) fflush( log );
}
#endif
double deltaIdeal = angle - idealAngle; double deltaIdeal = angle - idealAngle;
double deltaIdeal2 = deltaIdeal*deltaIdeal; double deltaIdeal2 = deltaIdeal*deltaIdeal;
double deltaIdeal3 = deltaIdeal*deltaIdeal2; double deltaIdeal3 = deltaIdeal*deltaIdeal2;
...@@ -121,7 +114,7 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou ...@@ -121,7 +114,7 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou
} }
static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions, AmoebaAngleForce& amoebaAngleForce, static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions, AmoebaAngleForce& amoebaAngleForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) { std::vector<Vec3>& forces, double* energy ) {
int particle1, particle2, particle3; int particle1, particle2, particle3;
double idealAngle; double idealAngle;
...@@ -133,14 +126,6 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions ...@@ -133,14 +126,6 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions
double penticK = amoebaAngleForce.getAmoebaGlobalAnglePentic(); double penticK = amoebaAngleForce.getAmoebaGlobalAnglePentic();
double sexticK = amoebaAngleForce.getAmoebaGlobalAngleSextic(); double sexticK = amoebaAngleForce.getAmoebaGlobalAngleSextic();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaAngleForce: bond %d [%d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, idealAngle, quadraticK, cubicK, quarticK, penticK, sexticK );
(void) fflush( log );
}
#endif
double deltaR[2][3]; double deltaR[2][3];
double r2_0 = 0.0; double r2_0 = 0.0;
double r2_1 = 0.0; double r2_1 = 0.0;
...@@ -163,17 +148,10 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions ...@@ -163,17 +148,10 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions
double dot = deltaR[0][0]*deltaR[1][0] + deltaR[0][1]*deltaR[1][1] + deltaR[0][2]*deltaR[1][2]; double dot = deltaR[0][0]*deltaR[1][0] + deltaR[0][1]*deltaR[1][1] + deltaR[0][2]*deltaR[1][2];
double cosine = dot/sqrt(r2_0*r2_1); double cosine = dot/sqrt(r2_0*r2_1);
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "dot=%10.3e r2_0=%10.3e r2_1=%10.3e\n", dot, r2_0, r2_1 );
(void) fflush( log );
}
#endif
double dEdR; double dEdR;
double energyTerm; double energyTerm;
getPrefactorsGivenAngleCosine( cosine, idealAngle, quadraticK, cubicK, getPrefactorsGivenAngleCosine( cosine, idealAngle, quadraticK, cubicK,
quarticK, penticK, sexticK, &dEdR, &energyTerm, log ); quarticK, penticK, sexticK, &dEdR, &energyTerm );
double termA = -dEdR/(r2_0*rp); double termA = -dEdR/(r2_0*rp);
double termC = dEdR/(r2_1*rp); double termC = dEdR/(r2_1*rp);
...@@ -203,7 +181,7 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions ...@@ -203,7 +181,7 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions
} }
static void computeAmoebaAngleForces( Context& context, AmoebaAngleForce& amoebaAngleForce, static void computeAmoebaAngleForces( Context& context, AmoebaAngleForce& amoebaAngleForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy ) {
// get positions and zero forces // get positions and zero forces
...@@ -219,51 +197,30 @@ static void computeAmoebaAngleForces( Context& context, AmoebaAngleForce& amoeba ...@@ -219,51 +197,30 @@ static void computeAmoebaAngleForces( Context& context, AmoebaAngleForce& amoeba
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaAngleForce.getNumAngles(); ii++ ){ for( int ii = 0; ii < amoebaAngleForce.getNumAngles(); ii++ ){
computeAmoebaAngleForce(ii, positions, amoebaAngleForce, expectedForces, expectedEnergy, log ); computeAmoebaAngleForce(ii, positions, amoebaAngleForce, expectedForces, expectedEnergy );
} }
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaAngleForces: 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; return;
} }
void compareWithExpectedForceAndEnergy( Context& context, AmoebaAngleForce& amoebaAngleForce, void compareWithExpectedForceAndEnergy( Context& context, AmoebaAngleForce& amoebaAngleForce,
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;
computeAmoebaAngleForces( context, amoebaAngleForce, expectedForces, &expectedEnergy, log ); computeAmoebaAngleForces( context, amoebaAngleForce, 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, "computeAmoebaAngleForces: 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 = 3; int numberOfParticles = 3;
...@@ -298,7 +255,7 @@ void testOneAngle( FILE* log ) { ...@@ -298,7 +255,7 @@ void testOneAngle( FILE* log ) {
positions[2] = Vec3(0, 0, 1); positions[2] = Vec3(0, 0, 1);
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log ); compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle" );
// Try changing the angle parameters and make sure it's still correct. // Try changing the angle parameters and make sure it's still correct.
...@@ -306,14 +263,14 @@ void testOneAngle( FILE* log ) { ...@@ -306,14 +263,14 @@ void testOneAngle( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log ); compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle" );
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaAngleForce->updateParametersInContext(context); amoebaAngleForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log ); compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle" );
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -322,8 +279,7 @@ int main(int argc, char* argv[]) { ...@@ -322,8 +279,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testOneAngle();
testOneAngle( log );
} 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/cuda/tests/TestCudaAmoebaBondForce.cpp
View file @ a568bb12
...@@ -84,7 +84,7 @@ static void computeAmoebaBondForce(int bondIndex, std::vector<Vec3>& positions, ...@@ -84,7 +84,7 @@ 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
...@@ -102,47 +102,24 @@ static void computeAmoebaBondForces( Context& context, AmoebaBondForce& amoebaBo ...@@ -102,47 +102,24 @@ static void computeAmoebaBondForces( Context& context, AmoebaBondForce& amoebaBo
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;
...@@ -169,10 +146,10 @@ void testOneBond( FILE* log ) { ...@@ -169,10 +146,10 @@ void testOneBond( FILE* log ) {
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;
...@@ -203,7 +180,7 @@ void testTwoBond( FILE* log ) { ...@@ -203,7 +180,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,14 +189,14 @@ void testTwoBond( FILE* log ) { ...@@ -212,14 +189,14 @@ 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 argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -228,8 +205,7 @@ int main(int argc, char* argv[]) { ...@@ -228,8 +205,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testTwoBond();
testTwoBond( log );
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl; std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp
View file @ a568bb12
...@@ -273,7 +273,7 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce ...@@ -273,7 +273,7 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
system.addForce(amoebaGeneralizedKirkwoodForce); system.addForce(amoebaGeneralizedKirkwoodForce);
} }
   
static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& forces, double& energy, FILE* log) { static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& forces, double& energy) {
std::vector<Vec3> positions(context.getSystem().getNumParticles()); std::vector<Vec3> positions(context.getSystem().getNumParticles());
   
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03 ); positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03 );
...@@ -294,7 +294,7 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& ...@@ -294,7 +294,7 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>&
// setup for villin // setup for villin
   
static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::PolarizationType polarizationType, static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::PolarizationType polarizationType,
int includeCavityTerm, std::vector<Vec3>& forces, double& energy, FILE* log ){ int includeCavityTerm, std::vector<Vec3>& forces, double& energy ){
   
// beginning of Multipole setup // beginning of Multipole setup
   
...@@ -6978,45 +6978,7 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari ...@@ -6978,45 +6978,7 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
   
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 );
} }
...@@ -7027,47 +6989,7 @@ static void compareForcesEnergy( std::string& testName, double expectedEnergy, d ...@@ -7027,47 +6989,7 @@ static void compareForcesEnergy( std::string& testName, double expectedEnergy, d
   
static void compareForceNormsEnergy( std::string& testName, double expectedEnergy, double energy, static void compareForceNormsEnergy( std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces, std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) { const std::vector<Vec3>& forces, double tolerance ) {
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
if( log ){
double conversion = 1.0/4.184;
double energyAbsDiff = fabs( expectedEnergy - energy );
double energyRelDiff = 2.0*energyAbsDiff/( fabs( expectedEnergy ) + fabs( energy ) + 1.0e-08 );
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e absDiff=%15.7e relDiff=%15.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy,
conversion*energyAbsDiff, conversion*energyRelDiff );
if( conversion != 1.0 )conversion *= -0.1;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2] );
double norm = sqrt( forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2] );
double absDiff = fabs( (norm - expectedNorm) );
double relDiff = 2.0*absDiff/(fabs( norm ) + fabs( expectedNorm ) + 1.0e-08);
(void) fprintf( log, "%6u %15.7e %15.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e] %15.7e %15.7e\n", ii,
fabs(conversion)*absDiff, relDiff,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2],
fabs(conversion)*expectedNorm, fabs(conversion)*norm );
}
(void) fflush( log );
conversion = 1.0;
(void) fprintf( log, "\n%s: expected energy=%14.7e %14.7e no conversion\n", testName.c_str(), conversion*expectedEnergy, conversion*energy );
if( conversion != 1.0 )conversion = -1.0;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] + double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] + expectedForces[ii][1]*expectedForces[ii][1] +
...@@ -7096,7 +7018,7 @@ static void compareForceNormsEnergy( std::string& testName, double expectedEnerg ...@@ -7096,7 +7018,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";
   
...@@ -7109,7 +7031,7 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) { ...@@ -7109,7 +7031,7 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) {
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0); setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
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;
...@@ -7124,12 +7046,12 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) { ...@@ -7124,12 +7046,12 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) {
expectedForces[7] = Vec3( 9.7274235e+00, -6.3130458e+01, 1.4949024e+02 ); expectedForces[7] = Vec3( 9.7274235e+00, -6.3130458e+01, 1.4949024e+02 );
   
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
   
// test GK mutual polarization for system comprised of two ammonia molecules // test GK mutual polarization for system comprised of two ammonia molecules
   
static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) { static void testGeneralizedKirkwoodAmmoniaMutualPolarization() {
   
std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarization"; std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarization";
   
...@@ -7142,7 +7064,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) { ...@@ -7142,7 +7064,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) {
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 0); setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
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;
...@@ -7157,13 +7079,13 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) { ...@@ -7157,13 +7079,13 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) {
expectedForces[7] = Vec3( 5.3895456e+00, -7.7131137e+01, 1.5826273e+02 ); expectedForces[7] = Vec3( 5.3895456e+00, -7.7131137e+01, 1.5826273e+02 );
   
double tolerance = 2.0e-04; double tolerance = 2.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";
   
...@@ -7176,7 +7098,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE ...@@ -7176,7 +7098,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 1); setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 1);
LangevinIntegrator integrator(0.0, 0.1, 0.01); LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA")); Context context(system, integrator, Platform::getPlatformByName("CUDA"));
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;
...@@ -7191,7 +7113,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE ...@@ -7191,7 +7113,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE
expectedForces[7] = Vec3( 7.9205382e+00, -7.3716473e+01, 1.5960993e+02 ); expectedForces[7] = Vec3( 7.9205382e+00, -7.3716473e+01, 1.5960993e+02 );
   
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance);
// Try changing the particle parameters and make sure it's still correct. // Try changing the particle parameters and make sure it's still correct.
...@@ -7208,7 +7130,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE ...@@ -7208,7 +7130,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;
...@@ -7216,12 +7138,12 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE ...@@ -7216,12 +7138,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";
   
...@@ -7229,7 +7151,7 @@ static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) { ...@@ -7229,7 +7151,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;
...@@ -7840,12 +7762,12 @@ static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) { ...@@ -7840,12 +7762,12 @@ static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) {
} }
   
double tolerance = 1.0e-05; double tolerance = 1.0e-05;
compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance);
} }
   
// test GK mutual polarization for villin system // test GK mutual polarization for villin system
   
static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) { static void testGeneralizedKirkwoodVillinMutualPolarization( ) {
   
std::string testName = "testGeneralizedKirkwoodVillinMutualPolarization"; std::string testName = "testGeneralizedKirkwoodVillinMutualPolarization";
   
...@@ -7853,7 +7775,7 @@ static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) { ...@@ -7853,7 +7775,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;
...@@ -8464,7 +8386,7 @@ static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) { ...@@ -8464,7 +8386,7 @@ static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) {
} }
   
double tolerance = 1.0e-05; double tolerance = 1.0e-05;
compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
   
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -8474,17 +8396,15 @@ int main(int argc, char* argv[]) { ...@@ -8474,17 +8396,15 @@ int main(int argc, char* argv[]) {
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
   
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
   
testGeneralizedKirkwoodAmmoniaDirectPolarization( log ); testGeneralizedKirkwoodAmmoniaDirectPolarization();
testGeneralizedKirkwoodAmmoniaMutualPolarization( log ); testGeneralizedKirkwoodAmmoniaMutualPolarization();
testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( log ); testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm();
   
testGeneralizedKirkwoodVillinDirectPolarization( log ); testGeneralizedKirkwoodVillinDirectPolarization();
testGeneralizedKirkwoodVillinMutualPolarization( log ); testGeneralizedKirkwoodVillinMutualPolarization();
   
} 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/cuda/tests/TestCudaAmoebaInPlaneAngleForce.cpp
View file @ a568bb12
...@@ -78,7 +78,7 @@ static double dotVector3( double* vectorX, double* vectorY ){ ...@@ -78,7 +78,7 @@ static double dotVector3( double* vectorX, double* vectorY ){
static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInPlaneAngle, double quadraticK, double cubicK, static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInPlaneAngle, double quadraticK, double cubicK,
double quarticK, double penticK, double sexticK, double quarticK, double penticK, double sexticK,
double* dEdR, double* energyTerm, FILE* log ) { double* dEdR, double* energyTerm ) {
double angle; double angle;
if( cosine >= 1.0 ){ if( cosine >= 1.0 ){
...@@ -89,13 +89,6 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP ...@@ -89,13 +89,6 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP
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;
double deltaIdeal3 = deltaIdeal*deltaIdeal2; double deltaIdeal3 = deltaIdeal*deltaIdeal2;
...@@ -122,7 +115,7 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP ...@@ -122,7 +115,7 @@ 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;
...@@ -134,14 +127,6 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -134,14 +127,6 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
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
// AP = A - P // AP = A - P
...@@ -181,14 +166,6 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -181,14 +166,6 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
double rAp2 = dotVector3( deltaR[AP], deltaR[AP] ); double rAp2 = dotVector3( deltaR[AP], deltaR[AP] );
double rCp2 = dotVector3( deltaR[CP], deltaR[CP] ); double rCp2 = dotVector3( deltaR[CP], deltaR[CP] );
if( rAp2 <= 0.0 && rCp2 <= 0.0 ){ if( rAp2 <= 0.0 && rCp2 <= 0.0 ){
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaInPlaneAngleForce: rAp2 or rCp2 <= 0.0\n" );
(void) fflush( log );
}
#endif
return;
} }
crossProductVector3( deltaR[CP], deltaR[AP], deltaR[M] ); crossProductVector3( deltaR[CP], deltaR[AP], deltaR[M] );
...@@ -205,7 +182,7 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -205,7 +182,7 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
double dEdR; double dEdR;
double energyTerm; double energyTerm;
getPrefactorsGivenInPlaneAngleCosine( cosine, idealInPlaneAngle, quadraticK, cubicK, getPrefactorsGivenInPlaneAngleCosine( cosine, idealInPlaneAngle, quadraticK, cubicK,
quarticK, penticK, sexticK, &dEdR, &energyTerm, log ); quarticK, penticK, sexticK, &dEdR, &energyTerm );
double termA = -dEdR/(rAp2*rm); double termA = -dEdR/(rAp2*rm);
double termC = dEdR/(rCp2*rm); double termC = dEdR/(rCp2*rm);
...@@ -280,7 +257,7 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po ...@@ -280,7 +257,7 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
} }
static void computeAmoebaInPlaneAngleForces( Context& context, AmoebaInPlaneAngleForce& amoebaInPlaneAngleForce, static void computeAmoebaInPlaneAngleForces( Context& context, AmoebaInPlaneAngleForce& amoebaInPlaneAngleForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy ) {
// get positions and zero forces // get positions and zero forces
...@@ -296,50 +273,26 @@ static void computeAmoebaInPlaneAngleForces( Context& context, AmoebaInPlaneAngl ...@@ -296,50 +273,26 @@ static void computeAmoebaInPlaneAngleForces( Context& context, AmoebaInPlaneAngl
*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;
...@@ -375,7 +328,7 @@ void testOneAngle( FILE* log ) { ...@@ -375,7 +328,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.
...@@ -383,14 +336,14 @@ void testOneAngle( FILE* log ) { ...@@ -383,14 +336,14 @@ void testOneAngle( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log ); compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle" );
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaInPlaneAngleForce->updateParametersInContext(context); amoebaInPlaneAngleForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log ); compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle" );
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -399,8 +352,7 @@ int main(int argc, char* argv[]) { ...@@ -399,8 +352,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testOneAngle();
testOneAngle( NULL );
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl; std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaMultipoleForce.cpp
View file @ a568bb12
...@@ -294,45 +294,7 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& ...@@ -294,45 +294,7 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>&
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 );
} }
...@@ -343,46 +305,7 @@ static void compareForcesEnergy( std::string& testName, double expectedEnergy, d ...@@ -343,46 +305,7 @@ static void compareForcesEnergy( std::string& testName, double expectedEnergy, d
static void compareForceNormsEnergy( std::string& testName, double expectedEnergy, double energy, static void compareForceNormsEnergy( std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces, std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) { const std::vector<Vec3>& forces, double tolerance ) {
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
if( log ){
double conversion = 1.0/4.184;
double energyAbsDiff = fabs( expectedEnergy - energy );
double energyRelDiff = 2.0*energyAbsDiff/( fabs( expectedEnergy ) + fabs( energy ) + 1.0e-08 );
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e absDiff=%15.7e relDiff=%15.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy,
conversion*energyAbsDiff, conversion*energyRelDiff );
if( conversion != 1.0 )conversion *= -0.1;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2] );
double norm = sqrt( forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2] );
double absDiff = fabs( (norm - expectedNorm) );
double relDiff = 2.0*absDiff/(fabs( norm ) + fabs( expectedNorm ) + 1.0e-08);
(void) fprintf( log, "%6u %15.7e %15.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e] %15.7e %15.7e\n", ii,
fabs(conversion)*absDiff, relDiff,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2],
fabs(conversion)*expectedNorm, fabs(conversion)*norm );
}
(void) fflush( log );
conversion = 1.0;
(void) fprintf( log, "\n%s: expected energy=%14.7e %14.7e no conversion\n", testName.c_str(), conversion*expectedEnergy, conversion*energy );
if( conversion != 1.0 )conversion = -1.0;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] + double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] + expectedForces[ii][1]*expectedForces[ii][1] +
...@@ -411,7 +334,7 @@ static void compareForceNormsEnergy( std::string& testName, double expectedEnerg ...@@ -411,7 +334,7 @@ static void compareForceNormsEnergy( std::string& testName, double expectedEnerg
// test multipole direct polarization for system comprised of two ammonia molecules; no cutoff // test multipole direct polarization for system comprised of two ammonia molecules; no cutoff
static void testMultipoleAmmoniaDirectPolarization( FILE* log ) { static void testMultipoleAmmoniaDirectPolarization() {
std::string testName = "testMultipoleAmmoniaDirectPolarization"; std::string testName = "testMultipoleAmmoniaDirectPolarization";
...@@ -442,12 +365,12 @@ static void testMultipoleAmmoniaDirectPolarization( FILE* log ) { ...@@ -442,12 +365,12 @@ static void testMultipoleAmmoniaDirectPolarization( FILE* log ) {
expectedForces[7] = Vec3( 4.1453480e+01, 1.6842405e+01, 1.6409513e+00 ); expectedForces[7] = Vec3( 4.1453480e+01, 1.6842405e+01, 1.6409513e+00 );
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
// test multipole mutual polarization for system comprised of two ammonia molecules; no cutoff // test multipole mutual polarization for system comprised of two ammonia molecules; no cutoff
static void testMultipoleAmmoniaMutualPolarization( FILE* log ) { static void testMultipoleAmmoniaMutualPolarization( ) {
std::string testName = "testMultipoleAmmoniaMutualPolarization"; std::string testName = "testMultipoleAmmoniaMutualPolarization";
...@@ -478,7 +401,7 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) { ...@@ -478,7 +401,7 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) {
expectedForces[7] = Vec3( 4.2293601e+01, 1.7186738e+01, 1.3017270e+00 ); expectedForces[7] = Vec3( 4.2293601e+01, 1.7186738e+01, 1.3017270e+00 );
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
// Try changing the particle parameters and make sure it's still correct. // Try changing the particle parameters and make sure it's still correct.
...@@ -500,7 +423,7 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) { ...@@ -500,7 +423,7 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) {
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 );
for (int i = 0; i < numberOfParticles; i++) for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], tolerance); ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], tolerance);
} }
...@@ -510,7 +433,7 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) { ...@@ -510,7 +433,7 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) {
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaMultipoleForce->updateParametersInContext(context); amoebaMultipoleForce->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 );
} }
// setup for box of 4 water molecules -- used to test PME // setup for box of 4 water molecules -- used to test PME
...@@ -518,7 +441,7 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) { ...@@ -518,7 +441,7 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) {
static void setupAndGetForcesEnergyMultipoleWater( AmoebaMultipoleForce::NonbondedMethod nonbondedMethod, static void setupAndGetForcesEnergyMultipoleWater( AmoebaMultipoleForce::NonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::PolarizationType polarizationType, AmoebaMultipoleForce::PolarizationType polarizationType,
double cutoff, int inputPmeGridDimension, std::vector<Vec3>& forces, double cutoff, int inputPmeGridDimension, std::vector<Vec3>& forces,
double& energy, FILE* log ){ double& energy){
// beginning of Multipole setup // beginning of Multipole setup
...@@ -671,7 +594,7 @@ static void setupAndGetForcesEnergyMultipoleWater( AmoebaMultipoleForce::Nonbond ...@@ -671,7 +594,7 @@ static void setupAndGetForcesEnergyMultipoleWater( AmoebaMultipoleForce::Nonbond
// test multipole direct polarization using PME for box of water // test multipole direct polarization using PME for box of water
static void testMultipoleWaterPMEDirectPolarization( FILE* log ) { static void testMultipoleWaterPMEDirectPolarization( ) {
std::string testName = "testMultipoleWaterDirectPolarization"; std::string testName = "testMultipoleWaterDirectPolarization";
...@@ -682,7 +605,7 @@ static void testMultipoleWaterPMEDirectPolarization( FILE* log ) { ...@@ -682,7 +605,7 @@ static void testMultipoleWaterPMEDirectPolarization( FILE* log ) {
double energy; double energy;
setupAndGetForcesEnergyMultipoleWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Direct, setupAndGetForcesEnergyMultipoleWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Direct,
cutoff, inputPmeGridDimension, forces, energy, log ); cutoff, inputPmeGridDimension, forces, energy );
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 6.4585115e-01; double expectedEnergy = 6.4585115e-01;
...@@ -701,12 +624,12 @@ static void testMultipoleWaterPMEDirectPolarization( FILE* log ) { ...@@ -701,12 +624,12 @@ static void testMultipoleWaterPMEDirectPolarization( FILE* log ) {
expectedForces[11] = Vec3( 1.2523841e+00, -1.9794292e+00, -3.4670129e+00 ); expectedForces[11] = Vec3( 1.2523841e+00, -1.9794292e+00, -3.4670129e+00 );
double tolerance = 1.0e-03; double tolerance = 1.0e-03;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
// test multipole mutual polarization using PME for box of water // test multipole mutual polarization using PME for box of water
static void testMultipoleWaterPMEMutualPolarization( FILE* log ) { static void testMultipoleWaterPMEMutualPolarization( ) {
std::string testName = "testMultipoleWaterMutualPolarization"; std::string testName = "testMultipoleWaterMutualPolarization";
...@@ -717,7 +640,7 @@ static void testMultipoleWaterPMEMutualPolarization( FILE* log ) { ...@@ -717,7 +640,7 @@ static void testMultipoleWaterPMEMutualPolarization( FILE* log ) {
double energy; double energy;
setupAndGetForcesEnergyMultipoleWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual, setupAndGetForcesEnergyMultipoleWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual,
cutoff, inputPmeGridDimension, forces, energy, log ); cutoff, inputPmeGridDimension, forces, energy );
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 6.5029855e-01; double expectedEnergy = 6.5029855e-01;
...@@ -736,12 +659,12 @@ static void testMultipoleWaterPMEMutualPolarization( FILE* log ) { ...@@ -736,12 +659,12 @@ static void testMultipoleWaterPMEMutualPolarization( FILE* log ) {
expectedForces[11] = Vec3( 1.2467701e+00, -1.9832979e+00, -3.4684052e+00 ); expectedForces[11] = Vec3( 1.2467701e+00, -1.9832979e+00, -3.4684052e+00 );
double tolerance = 1.0e-03; double tolerance = 1.0e-03;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
// check validation of traceless/symmetric quadrupole tensor // check validation of traceless/symmetric quadrupole tensor
static void testQuadrupoleValidation( FILE* log ){ static void testQuadrupoleValidation( ){
std::string testName = "checkQuadrupoleValidation"; std::string testName = "checkQuadrupoleValidation";
...@@ -958,7 +881,7 @@ static void testQuadrupoleValidation( FILE* log ){ ...@@ -958,7 +881,7 @@ static void testQuadrupoleValidation( FILE* log ){
static void setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::NonbondedMethod nonbondedMethod, static void setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::NonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::PolarizationType polarizationType, AmoebaMultipoleForce::PolarizationType polarizationType,
double cutoff, int inputPmeGridDimension, std::string testName, double cutoff, int inputPmeGridDimension, std::string testName,
std::vector<Vec3>& forces, double& energy, FILE* log ){ std::vector<Vec3>& forces, double& energy ){
// beginning of Multipole setup // beginning of Multipole setup
...@@ -1149,7 +1072,7 @@ static void setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce:: ...@@ -1149,7 +1072,7 @@ static void setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::
// test multipole mutual polarization using PME for system comprised of 2 ions and 2 waters // test multipole mutual polarization using PME for system comprised of 2 ions and 2 waters
static void testMultipoleIonsAndWaterPMEDirectPolarization( FILE* log ) { static void testMultipoleIonsAndWaterPMEDirectPolarization( ) {
std::string testName = "testMultipoleIonsAndWaterDirectPolarization"; std::string testName = "testMultipoleIonsAndWaterDirectPolarization";
...@@ -1161,7 +1084,7 @@ static void testMultipoleIonsAndWaterPMEDirectPolarization( FILE* log ) { ...@@ -1161,7 +1084,7 @@ static void testMultipoleIonsAndWaterPMEDirectPolarization( FILE* log ) {
double energy; double energy;
setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Direct, setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Direct,
cutoff, inputPmeGridDimension, testName, forces, energy, log ); cutoff, inputPmeGridDimension, testName, forces, energy );
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
...@@ -1177,13 +1100,13 @@ static void testMultipoleIonsAndWaterPMEDirectPolarization( FILE* log ) { ...@@ -1177,13 +1100,13 @@ static void testMultipoleIonsAndWaterPMEDirectPolarization( FILE* log ) {
expectedForces[7] = Vec3( 4.0644209e+00, -3.3666305e+00, -1.7022384e+00 ); expectedForces[7] = Vec3( 4.0644209e+00, -3.3666305e+00, -1.7022384e+00 );
double tolerance = 5.0e-04; double tolerance = 5.0e-04;
compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
// test multipole mutual polarization using PME for system comprised of 2 ions and 2 waters // test multipole mutual polarization using PME for system comprised of 2 ions and 2 waters
static void testMultipoleIonsAndWaterPMEMutualPolarization( FILE* log ) { static void testMultipoleIonsAndWaterPMEMutualPolarization( ) {
std::string testName = "testMultipoleIonsAndWaterMutualPolarization"; std::string testName = "testMultipoleIonsAndWaterMutualPolarization";
...@@ -1197,7 +1120,7 @@ static void testMultipoleIonsAndWaterPMEMutualPolarization( FILE* log ) { ...@@ -1197,7 +1120,7 @@ static void testMultipoleIonsAndWaterPMEMutualPolarization( FILE* log ) {
std::vector<Vec3> inputGrid; std::vector<Vec3> inputGrid;
setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual, setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual,
cutoff, inputPmeGridDimension, testName, forces, energy, log ); cutoff, inputPmeGridDimension, testName, forces, energy );
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
...@@ -1213,9 +1136,9 @@ static void testMultipoleIonsAndWaterPMEMutualPolarization( FILE* log ) { ...@@ -1213,9 +1136,9 @@ static void testMultipoleIonsAndWaterPMEMutualPolarization( FILE* log ) {
expectedForces[7] = Vec3( 4.0622614e+00, -3.3687594e+00, -1.6986575e+00 ); expectedForces[7] = Vec3( 4.0622614e+00, -3.3687594e+00, -1.6986575e+00 );
//double tolerance = 1.0e-03; //double tolerance = 1.0e-03;
//compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); //compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
double tolerance = 5.0e-04; double tolerance = 5.0e-04;
compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
...@@ -1227,7 +1150,7 @@ static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::No ...@@ -1227,7 +1150,7 @@ static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::No
std::vector<Vec3>& forces, double& energy, std::vector<Vec3>& forces, double& energy,
std::vector< double >& outputMultipoleMoments, std::vector< double >& outputMultipoleMoments,
std::vector< Vec3 >& inputGrid, std::vector< Vec3 >& inputGrid,
std::vector< double >& outputGridPotential, FILE* log ){ std::vector< double >& outputGridPotential ){
// beginning of Multipole setup // beginning of Multipole setup
...@@ -2025,7 +1948,7 @@ static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::No ...@@ -2025,7 +1948,7 @@ static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::No
// test multipole mutual polarization using PME for box of water // test multipole mutual polarization using PME for box of water
static void testPMEMutualPolarizationLargeWater( FILE* log ) { static void testPMEMutualPolarizationLargeWater( ) {
std::string testName = "testPMEMutualPolarizationLargeWater"; std::string testName = "testPMEMutualPolarizationLargeWater";
...@@ -2040,7 +1963,7 @@ static void testPMEMutualPolarizationLargeWater( FILE* log ) { ...@@ -2040,7 +1963,7 @@ static void testPMEMutualPolarizationLargeWater( FILE* log ) {
setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual, setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual,
cutoff, inputPmeGridDimension, testName, cutoff, inputPmeGridDimension, testName,
forces, energy, outputMultipoleMoments, inputGrid, outputGridPotential, log ); forces, energy, outputMultipoleMoments, inputGrid, outputGridPotential );
static std::vector<Vec3> expectedForces; // Static to work around bug in Visual Studio that makes compilation very very slow. static std::vector<Vec3> expectedForces; // Static to work around bug in Visual Studio that makes compilation very very slow.
expectedForces.resize(numberOfParticles); expectedForces.resize(numberOfParticles);
...@@ -2696,7 +2619,7 @@ static void testPMEMutualPolarizationLargeWater( FILE* log ) { ...@@ -2696,7 +2619,7 @@ static void testPMEMutualPolarizationLargeWater( FILE* log ) {
expectedForces[647] = Vec3( -4.7388806e+02, -5.5561844e+02, -8.5019295e+02 ); expectedForces[647] = Vec3( -4.7388806e+02, -5.5561844e+02, -8.5019295e+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 );
} }
...@@ -2736,7 +2659,7 @@ static void testParticleInducedDipoles() { ...@@ -2736,7 +2659,7 @@ static void testParticleInducedDipoles() {
// test computation of system multipole moments // test computation of system multipole moments
static void testSystemMultipoleMoments( FILE* log ) { static void testSystemMultipoleMoments( ) {
std::string testName = "testSystemMultipoleMoments"; std::string testName = "testSystemMultipoleMoments";
...@@ -2752,7 +2675,7 @@ static void testSystemMultipoleMoments( FILE* log ) { ...@@ -2752,7 +2675,7 @@ static void testSystemMultipoleMoments( FILE* log ) {
setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual, setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual,
cutoff, inputPmeGridDimension, testName, cutoff, inputPmeGridDimension, testName,
forces, energy, outputMultipoleMoments, inputGrid, outputGridPotential, log ); forces, energy, outputMultipoleMoments, inputGrid, outputGridPotential );
std::vector<double> tinkerMoments(4); std::vector<double> tinkerMoments(4);
...@@ -2771,9 +2694,6 @@ static void testSystemMultipoleMoments( FILE* log ) { ...@@ -2771,9 +2694,6 @@ static void testSystemMultipoleMoments( FILE* log ) {
// tinkerMoments[12] = 4.3292490e-02; // tinkerMoments[12] = 4.3292490e-02;
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
if( log ){
(void) fprintf( log, "%s RelativeDifference Tinker OpenMM\n", testName.c_str() );
}
for( unsigned int ii = 0; ii < tinkerMoments.size(); ii++ ){ for( unsigned int ii = 0; ii < tinkerMoments.size(); ii++ ){
double difference; double difference;
if( fabs( tinkerMoments[ii] ) > 0.0 ){ if( fabs( tinkerMoments[ii] ) > 0.0 ){
...@@ -2781,10 +2701,6 @@ static void testSystemMultipoleMoments( FILE* log ) { ...@@ -2781,10 +2701,6 @@ static void testSystemMultipoleMoments( FILE* log ) {
} else { } else {
difference = fabs( outputMultipoleMoments[ii] - tinkerMoments[ii] ); difference = fabs( outputMultipoleMoments[ii] - tinkerMoments[ii] );
} }
if( log ){
(void) fprintf( log, "%2d %15.7e %15.7e %15.7e\n", ii, difference, tinkerMoments[ii], outputMultipoleMoments[ii] );
}
if( difference > tolerance ){ if( difference > tolerance ){
std::stringstream details; std::stringstream details;
details << testName << "Multipole moment " << ii << " does not agree w/ TINKER computed moments: OpenMM=" << outputMultipoleMoments[ii]; details << testName << "Multipole moment " << ii << " does not agree w/ TINKER computed moments: OpenMM=" << outputMultipoleMoments[ii];
...@@ -2798,7 +2714,7 @@ static void testSystemMultipoleMoments( FILE* log ) { ...@@ -2798,7 +2714,7 @@ static void testSystemMultipoleMoments( FILE* log ) {
// test computation of multipole potential on a grid // test computation of multipole potential on a grid
static void testMultipoleGridPotential( FILE* log ) { static void testMultipoleGridPotential( ) {
std::string testName = "testMultipoleGridPotential"; std::string testName = "testMultipoleGridPotential";
...@@ -2847,7 +2763,7 @@ static void testMultipoleGridPotential( FILE* log ) { ...@@ -2847,7 +2763,7 @@ static void testMultipoleGridPotential( FILE* log ) {
setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual, setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::PME, AmoebaMultipoleForce::Mutual,
cutoff, inputPmeGridDimension, testName, forces, energy, cutoff, inputPmeGridDimension, testName, forces, energy,
outputMultipoleMoments, inputGrid, outputGridPotential, log ); outputMultipoleMoments, inputGrid, outputGridPotential );
// TINKER computed grid values // TINKER computed grid values
...@@ -2943,7 +2859,7 @@ static void testNoQuadrupoles(bool usePme) { ...@@ -2943,7 +2859,7 @@ static void testNoQuadrupoles(bool usePme) {
double energy3; double energy3;
getForcesEnergyMultipoleAmmonia(context2, forces3, energy3); getForcesEnergyMultipoleAmmonia(context2, forces3, energy3);
double tolerance = 1e-5; double tolerance = 1e-5;
compareForcesEnergy(testName, energy2, energy3, forces2, forces3, tolerance, NULL); compareForcesEnergy(testName, energy2, energy3, forces2, forces3, tolerance);
// If we try to set a quadrupole, that should produce and exception. // If we try to set a quadrupole, that should produce and exception.
...@@ -3144,31 +3060,29 @@ int main(int argc, char* argv[]) { ...@@ -3144,31 +3060,29 @@ int main(int argc, char* argv[]) {
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL;
// tests using two ammonia molecules // tests using two ammonia molecules
// test direct polarization, no cutoff // test direct polarization, no cutoff
testMultipoleAmmoniaDirectPolarization( log ); testMultipoleAmmoniaDirectPolarization();
// test mutual polarization, no cutoff // test mutual polarization, no cutoff
testMultipoleAmmoniaMutualPolarization( log ); testMultipoleAmmoniaMutualPolarization();
// test multipole direct & mutual polarization using PME // test multipole direct & mutual polarization using PME
testMultipoleWaterPMEDirectPolarization( log ); testMultipoleWaterPMEDirectPolarization();
testMultipoleWaterPMEMutualPolarization( log ); testMultipoleWaterPMEMutualPolarization();
// check validation of traceless/symmetric quadrupole tensor // check validation of traceless/symmetric quadrupole tensor
testQuadrupoleValidation( log ); testQuadrupoleValidation();
// system w/ 2 ions and 2 water molecules // system w/ 2 ions and 2 water molecules
testMultipoleIonsAndWaterPMEMutualPolarization( log ); testMultipoleIonsAndWaterPMEMutualPolarization();
testMultipoleIonsAndWaterPMEDirectPolarization( log ); testMultipoleIonsAndWaterPMEDirectPolarization();
// test querying induced dipoles // test querying induced dipoles
...@@ -3176,14 +3090,14 @@ int main(int argc, char* argv[]) { ...@@ -3176,14 +3090,14 @@ int main(int argc, char* argv[]) {
// test computation of system multipole moments // test computation of system multipole moments
testSystemMultipoleMoments( log ); testSystemMultipoleMoments();
// test computation of grid potential // test computation of grid potential
testMultipoleGridPotential( log ); testMultipoleGridPotential();
// large box of water // large box of water
testPMEMutualPolarizationLargeWater( log ); testPMEMutualPolarizationLargeWater();
testNoQuadrupoles(false); testNoQuadrupoles(false);
testNoQuadrupoles(true); testNoQuadrupoles(true);
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaOutOfPlaneBendForce.cpp
View file @ a568bb12
...@@ -78,7 +78,7 @@ static double dotVector3( double* vectorX, double* vectorY ){ ...@@ -78,7 +78,7 @@ static double dotVector3( double* vectorX, double* vectorY ){
static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce, static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) { std::vector<Vec3>& forces, double* energy ) {
double kAngleCubic = amoebaOutOfPlaneBendForce.getAmoebaGlobalOutOfPlaneBendCubic(); double kAngleCubic = amoebaOutOfPlaneBendForce.getAmoebaGlobalOutOfPlaneBendCubic();
...@@ -90,14 +90,6 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -90,14 +90,6 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
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 };
...@@ -139,14 +131,6 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -139,14 +131,6 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
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
double dt = angle; double dt = angle;
...@@ -242,7 +226,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -242,7 +226,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
} }
static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce, static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy ) {
// get positions and zero forces // get positions and zero forces
...@@ -258,51 +242,26 @@ static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlan ...@@ -258,51 +242,26 @@ static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlan
*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;
...@@ -334,7 +293,7 @@ void testOneOutOfPlaneBend( FILE* log ) { ...@@ -334,7 +293,7 @@ void testOneOutOfPlaneBend( FILE* log ) {
positions[3] = Vec3( 0.245568230E+02, 0.250215290E+02, 0.796852800E+01 ); positions[3] = Vec3( 0.245568230E+02, 0.250215290E+02, 0.796852800E+01 );
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log ); compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
// Try changing the bend parameters and make sure it's still correct. // Try changing the bend parameters and make sure it's still correct.
...@@ -342,17 +301,17 @@ void testOneOutOfPlaneBend( FILE* log ) { ...@@ -342,17 +301,17 @@ void testOneOutOfPlaneBend( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log ); compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaOutOfPlaneBendForce->updateParametersInContext(context); amoebaOutOfPlaneBendForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log ); compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
} }
void testOneOutOfPlaneBend2( FILE* log, int setId ) { void testOneOutOfPlaneBend2(int setId ) {
System system; System system;
int numberOfParticles = 4; int numberOfParticles = 4;
...@@ -436,21 +395,10 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -436,21 +395,10 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
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);
...@@ -459,12 +407,6 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -459,12 +407,6 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for( unsigned int ii = 0; ii < 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() );
...@@ -473,7 +415,7 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -473,7 +415,7 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
} }
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log ); compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend");
static int iter = 0; static int iter = 0;
static std::map<int,Vec3> totalForces; static std::map<int,Vec3> totalForces;
...@@ -496,7 +438,7 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -496,7 +438,7 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
forces.resize( numberOfParticles ); forces.resize( numberOfParticles );
double energy; double energy;
computeAmoebaOutOfPlaneBendForce( 0, positions, *amoebaOutOfPlaneBendForce, forces, &energy, log ); computeAmoebaOutOfPlaneBendForce( 0, positions, *amoebaOutOfPlaneBendForce, forces, &energy);
totalEnergy += energy; totalEnergy += energy;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){ for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){
...@@ -504,22 +446,6 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -504,22 +446,6 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
totalForces[particleIndices[ii]][jj] += forces[ii][jj]; totalForces[particleIndices[ii]][jj] += forces[ii][jj];
} }
} }
if( iter == 6 ){
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: energy=%14.7e\n", totalEnergy);
for( std::map<int,Vec3>::iterator ii = totalForces.begin(); ii != totalForces.end(); ii++ ){
int particleIndex = ii->first;
Vec3 forces = ii->second;
(void) fprintf( log, "%6d [%14.7e %14.7e %14.7e] \n", particleIndex,
forces[0], forces[1], forces[2] );
}
(void) fflush( log );
}
#endif
}
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -528,9 +454,8 @@ int main(int argc, char* argv[]) { ...@@ -528,9 +454,8 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL;
testOneOutOfPlaneBend( log ); testOneOutOfPlaneBend();
} 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/cuda/tests/TestCudaAmoebaPiTorsionForce.cpp
View file @ a568bb12
...@@ -78,21 +78,13 @@ static double dotVector3( double* vectorX, double* vectorY ){ ...@@ -78,21 +78,13 @@ static double dotVector3( double* vectorX, double* vectorY ){
static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& positions, AmoebaPiTorsionForce& amoebaPiTorsionForce, static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& positions, AmoebaPiTorsionForce& amoebaPiTorsionForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) { std::vector<Vec3>& forces, double* energy) {
int particle1, particle2, particle3, particle4, particle5, particle6; int particle1, particle2, particle3, particle4, particle5, particle6;
double kTorsion; double kTorsion;
amoebaPiTorsionForce.getPiTorsionParameters(bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion); amoebaPiTorsionForce.getPiTorsionParameters(bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion);
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForce: bond %d [%d %d %d %d %d %d] k=%10.3e\n",
bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion );
(void) fflush( log );
}
#endif
enum { AD, BD, EC, FC, P, Q, CP, DC, QD, T, U, TU, DP, QC, dT, dU, dP, dQ, dC1, dC2, dD1, dD2, LastDeltaIndex }; enum { AD, BD, EC, FC, P, Q, CP, DC, QD, T, U, TU, DP, QC, dT, dU, dP, dQ, dC1, dC2, dD1, dD2, LastDeltaIndex };
double deltaR[LastDeltaIndex][3]; double deltaR[LastDeltaIndex][3];
...@@ -213,7 +205,7 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit ...@@ -213,7 +205,7 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
} }
static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce, static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy ) {
// get positions and zero forces // get positions and zero forces
...@@ -229,50 +221,27 @@ static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce ...@@ -229,50 +221,27 @@ static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaPiTorsionForce.getNumPiTorsions(); ii++ ){ for( int ii = 0; ii < amoebaPiTorsionForce.getNumPiTorsions(); ii++ ){
computeAmoebaPiTorsionForce(ii, positions, amoebaPiTorsionForce, expectedForces, expectedEnergy, log ); computeAmoebaPiTorsionForce(ii, positions, amoebaPiTorsionForce, expectedForces, expectedEnergy);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
}
(void) fflush( log );
} }
#endif
return;
} }
void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce, void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
double tolerance, const std::string& idString, FILE* log) { double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaPiTorsionForces( context, amoebaPiTorsionForce, expectedForces, &expectedEnergy, log ); computeAmoebaPiTorsionForces( context, amoebaPiTorsionForce, expectedForces, &expectedEnergy );
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance ); ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
} }
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;
...@@ -301,7 +270,7 @@ void testOnePiTorsion( FILE* log ) { ...@@ -301,7 +270,7 @@ void testOnePiTorsion( FILE* log ) {
positions[5] = Vec3( 0.254124630E+02, 0.234691880E+02, 0.773335400E+01 ); positions[5] = Vec3( 0.254124630E+02, 0.234691880E+02, 0.773335400E+01 );
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log ); compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion");
// Try changing the torsion parameters and make sure it's still correct. // Try changing the torsion parameters and make sure it's still correct.
...@@ -309,14 +278,14 @@ void testOnePiTorsion( FILE* log ) { ...@@ -309,14 +278,14 @@ void testOnePiTorsion( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log ); compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion");
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaPiTorsionForce->updateParametersInContext(context); amoebaPiTorsionForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log ); compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion" );
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -325,8 +294,7 @@ int main(int argc, char* argv[]) { ...@@ -325,8 +294,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testOnePiTorsion();
testOnePiTorsion( log );
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl; std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaStretchBendForce.cpp
View file @ a568bb12
...@@ -79,21 +79,13 @@ static double dotVector3( double* vectorX, double* vectorY ){ ...@@ -79,21 +79,13 @@ static double dotVector3( double* vectorX, double* vectorY ){
static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaStretchBendForce& amoebaStretchBendForce, static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaStretchBendForce& amoebaStretchBendForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) { std::vector<Vec3>& forces, double* energy) {
int particle1, particle2, particle3; int particle1, particle2, particle3;
double abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend; double abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend;
amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend); amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend);
angleStretchBend *= RADIAN; angleStretchBend *= RADIAN;
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForce: bond %d [%d %d %d] ab=%10.3e cb=%10.3e angle=%10.3e k1=%10.3e k2=%10.3e\n",
bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend );
(void) fflush( log );
}
#endif
enum { A, B, C, LastAtomIndex }; enum { A, B, C, LastAtomIndex };
enum { AB, CB, CBxAB, ABxP, CBxP, LastDeltaIndex }; enum { AB, CB, CBxAB, ABxP, CBxP, LastDeltaIndex };
...@@ -188,18 +180,10 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -188,18 +180,10 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
forces[particle3][2] -= subForce[2][2]; forces[particle3][2] -= subForce[2][2];
*energy += dt*drkk; *energy += dt*drkk;
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForce: angle=%10.3e dt=%10.3e dr=%10.3e\n", angle, dt, dr );
(void) fflush( log );
}
#endif
return;
} }
static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendForce& amoebaStretchBendForce, static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) { std::vector<Vec3>& expectedForces, double* expectedEnergy) {
// get positions and zero forces // get positions and zero forces
...@@ -215,50 +199,26 @@ static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendF ...@@ -215,50 +199,26 @@ static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendF
*expectedEnergy = 0.0; *expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaStretchBendForce.getNumStretchBends(); ii++ ){ for( int ii = 0; ii < amoebaStretchBendForce.getNumStretchBends(); ii++ ){
computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy, log ); computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy );
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
}
(void) fflush( log );
} }
#endif
return;
} }
void compareWithExpectedForceAndEnergy( Context& context, AmoebaStretchBendForce& amoebaStretchBendForce, void compareWithExpectedForceAndEnergy( Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
double tolerance, const std::string& idString, FILE* log) { double tolerance, const std::string& idString) {
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
double expectedEnergy; double expectedEnergy;
computeAmoebaStretchBendForces( context, amoebaStretchBendForce, expectedForces, &expectedEnergy, log ); computeAmoebaStretchBendForces( context, amoebaStretchBendForce, expectedForces, &expectedEnergy );
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces(); const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance ); ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
} }
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance ); ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
} }
void testOneStretchBend( FILE* log ) { void testOneStretchBend() {
System system; System system;
int numberOfParticles = 3; int numberOfParticles = 3;
...@@ -288,7 +248,7 @@ void testOneStretchBend( FILE* log ) { ...@@ -288,7 +248,7 @@ void testOneStretchBend( FILE* log ) {
positions[2] = Vec3( 0.269573220E+02, 0.236108860E+02, 0.216376800E+01 ); positions[2] = Vec3( 0.269573220E+02, 0.236108860E+02, 0.216376800E+01 );
context.setPositions(positions); context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log ); compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend" );
// Try changing the stretch-bend parameters and make sure it's still correct. // Try changing the stretch-bend parameters and make sure it's still correct.
...@@ -296,14 +256,14 @@ void testOneStretchBend( FILE* log ) { ...@@ -296,14 +256,14 @@ void testOneStretchBend( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log ); compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend" );
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
} }
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaStretchBendForce->updateParametersInContext(context); amoebaStretchBendForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log ); compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend" );
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -312,8 +272,7 @@ int main(int argc, char* argv[]) { ...@@ -312,8 +272,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testOneStretchBend();
testOneStretchBend( log );
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl; std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaTorsionTorsionForce.cpp
View file @ a568bb12
...@@ -2584,7 +2584,7 @@ TorsionTorsionGrid& getTorsionGrid( int gridIndex ) { ...@@ -2584,7 +2584,7 @@ TorsionTorsionGrid& getTorsionGrid( int gridIndex ) {
return grids[gridIndex]; return grids[gridIndex];
} }
void testTorsionTorsion( FILE* log, int systemId ) { void testTorsionTorsion(int systemId) {
System system; System system;
int numberOfParticles = 6; int numberOfParticles = 6;
...@@ -2660,19 +2660,6 @@ void testTorsionTorsion( FILE* log, int systemId ) { ...@@ -2660,19 +2660,6 @@ void testTorsionTorsion( FILE* log, int systemId ) {
forces[ii][1] *= conversion; forces[ii][1] *= conversion;
forces[ii][2] *= conversion; forces[ii][2] *= conversion;
} }
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaTorsionTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2],
forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
double tolerance = 1.0e-03; double tolerance = 1.0e-03;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance ); ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
...@@ -2687,8 +2674,7 @@ int main(int argc, char* argv[]) { ...@@ -2687,8 +2674,7 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL; testTorsionTorsion(1);
testTorsionTorsion( log, 1 );
} catch(const std::exception& e) { } catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl; std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl; std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaVdwForce.cpp
View file @ a568bb12
...@@ -57,7 +57,7 @@ extern "C" void registerAmoebaCudaKernelFactories(); ...@@ -57,7 +57,7 @@ extern "C" void registerAmoebaCudaKernelFactories();
const double TOL = 1e-4; const double TOL = 1e-4;
void testVdw( FILE* log ) { void testVdw() {
System system; System system;
int numberOfParticles = 6; int numberOfParticles = 6;
...@@ -159,18 +159,6 @@ void testVdw( FILE* log ) { ...@@ -159,18 +159,6 @@ void testVdw( FILE* log ) {
forces[ii][2] *= conversion; forces[ii][2] *= conversion;
} }
expectedEnergy *= CalToJoule; expectedEnergy *= CalToJoule;
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaVdwForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
double tolerance = 1.0e-03; double tolerance = 1.0e-03;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance ); ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
...@@ -208,7 +196,7 @@ void testVdw( FILE* log ) { ...@@ -208,7 +196,7 @@ void testVdw( FILE* log ) {
} }
void setupAndGetForcesEnergyVdwAmmonia( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff, void setupAndGetForcesEnergyVdwAmmonia( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff,
double boxDimension, std::vector<Vec3>& forces, double& energy, FILE* log ){ double boxDimension, std::vector<Vec3>& forces, double& energy){
// beginning of Vdw setup // beginning of Vdw setup
...@@ -343,24 +331,7 @@ void setupAndGetForcesEnergyVdwAmmonia( const std::string& sigmaCombiningRule, c ...@@ -343,24 +331,7 @@ void setupAndGetForcesEnergyVdwAmmonia( const std::string& sigmaCombiningRule, c
void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy, void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces, std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance, FILE* log ) { std::vector<Vec3>& forces, double tolerance) {
#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
if( log ){
double conversion = 1.0/4.184;
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy );
conversion *= -0.1;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC_MOD( expectedForces[ii], forces[ii], tolerance, testName ); ASSERT_EQUAL_VEC_MOD( expectedForces[ii], forces[ii], tolerance, testName );
} }
...@@ -369,7 +340,7 @@ void compareForcesEnergy( std::string& testName, double expectedEnergy, double e ...@@ -369,7 +340,7 @@ void compareForcesEnergy( std::string& testName, double expectedEnergy, double e
// test VDW w/ sigmaRule=CubicMean and epsilonRule=HHG // test VDW w/ sigmaRule=CubicMean and epsilonRule=HHG
void testVdwAmmoniaCubicMeanHhg( FILE* log ) { void testVdwAmmoniaCubicMeanHhg() {
std::string testName = "testVdwAmmoniaCubicMeanHhg"; std::string testName = "testVdwAmmoniaCubicMeanHhg";
...@@ -379,7 +350,7 @@ void testVdwAmmoniaCubicMeanHhg( FILE* log ) { ...@@ -379,7 +350,7 @@ void testVdwAmmoniaCubicMeanHhg( FILE* log ) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy, log ); setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy );
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 4.8012258e+00; double expectedEnergy = 4.8012258e+00;
...@@ -394,12 +365,12 @@ void testVdwAmmoniaCubicMeanHhg( FILE* log ) { ...@@ -394,12 +365,12 @@ void testVdwAmmoniaCubicMeanHhg( FILE* log ) {
expectedForces[7] = Vec3( 1.6756544e+00, 3.2497316e-01, -1.7906832e-01 ); expectedForces[7] = Vec3( 1.6756544e+00, 3.2497316e-01, -1.7906832e-01 );
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
// test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic // test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
void testVdwAmmoniaArithmeticArithmetic( FILE* log ) { void testVdwAmmoniaArithmeticArithmetic() {
std::string testName = "testVdwAmmoniaArithmeticArithmetic"; std::string testName = "testVdwAmmoniaArithmeticArithmetic";
...@@ -409,7 +380,7 @@ void testVdwAmmoniaArithmeticArithmetic( FILE* log ) { ...@@ -409,7 +380,7 @@ void testVdwAmmoniaArithmeticArithmetic( FILE* log ) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia( "ARITHMETIC", "ARITHMETIC", cutoff, boxDimension, forces, energy, log ); setupAndGetForcesEnergyVdwAmmonia( "ARITHMETIC", "ARITHMETIC", cutoff, boxDimension, forces, energy );
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 4.2252403e+00; double expectedEnergy = 4.2252403e+00;
...@@ -424,12 +395,12 @@ void testVdwAmmoniaArithmeticArithmetic( FILE* log ) { ...@@ -424,12 +395,12 @@ void testVdwAmmoniaArithmeticArithmetic( FILE* log ) {
expectedForces[7] = Vec3( 2.3761408e+00, 4.6871961e-01, -2.4731607e-01 ); expectedForces[7] = Vec3( 2.3761408e+00, 4.6871961e-01, -2.4731607e-01 );
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
// test VDW w/ sigmaRule=Geometric and epsilonRule=Geometric // test VDW w/ sigmaRule=Geometric and epsilonRule=Geometric
void testVdwAmmoniaGeometricGeometric( FILE* log ) { void testVdwAmmoniaGeometricGeometric( ) {
std::string testName = "testVdwAmmoniaGeometricGeometric"; std::string testName = "testVdwAmmoniaGeometricGeometric";
...@@ -438,7 +409,7 @@ void testVdwAmmoniaGeometricGeometric( FILE* log ) { ...@@ -438,7 +409,7 @@ void testVdwAmmoniaGeometricGeometric( FILE* log ) {
double cutoff = 9000000.0; double cutoff = 9000000.0;
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia( "GEOMETRIC", "GEOMETRIC", cutoff, boxDimension, forces, energy, log ); setupAndGetForcesEnergyVdwAmmonia( "GEOMETRIC", "GEOMETRIC", cutoff, boxDimension, forces, energy );
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 2.5249914e+00; double expectedEnergy = 2.5249914e+00;
...@@ -453,10 +424,10 @@ void testVdwAmmoniaGeometricGeometric( FILE* log ) { ...@@ -453,10 +424,10 @@ void testVdwAmmoniaGeometricGeometric( FILE* log ) {
expectedForces[7] = Vec3( 1.8109211e+00, 3.5273117e-01, -1.9224723e-01 ); expectedForces[7] = Vec3( 1.8109211e+00, 3.5273117e-01, -1.9224723e-01 );
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
void testVdwAmmoniaCubicMeanHarmonic( FILE* log ) { void testVdwAmmoniaCubicMeanHarmonic( ) {
std::string testName = "testVdwAmmoniaCubicMeanHarmonic"; std::string testName = "testVdwAmmoniaCubicMeanHarmonic";
...@@ -465,7 +436,7 @@ void testVdwAmmoniaCubicMeanHarmonic( FILE* log ) { ...@@ -465,7 +436,7 @@ void testVdwAmmoniaCubicMeanHarmonic( FILE* log ) {
double cutoff = 9000000.0; double cutoff = 9000000.0;
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HARMONIC", cutoff, boxDimension, forces, energy, log ); setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HARMONIC", cutoff, boxDimension, forces, energy );
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 4.1369069e+00; double expectedEnergy = 4.1369069e+00;
...@@ -480,14 +451,14 @@ void testVdwAmmoniaCubicMeanHarmonic( FILE* log ) { ...@@ -480,14 +451,14 @@ void testVdwAmmoniaCubicMeanHarmonic( FILE* log ) {
expectedForces[7] = Vec3( 1.5080748e+00, 2.9058422e-01, -1.6274118e-01 ); expectedForces[7] = Vec3( 1.5080748e+00, 2.9058422e-01, -1.6274118e-01 );
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
// test w/ cutoff=0.25 nm; single ixn between two particles (0 and 6); force nonzero on // test w/ cutoff=0.25 nm; single ixn between two particles (0 and 6); force nonzero on
// particle 4 due to reduction applied to NH // particle 4 due to reduction applied to NH
// the distance between 0 and 6 is ~ 0.235 so the ixn is in the tapered region // the distance between 0 and 6 is ~ 0.235 so the ixn is in the tapered region
void testVdwTaper( FILE* log ) { void testVdwTaper( ) {
std::string testName = "testVdwTaper"; std::string testName = "testVdwTaper";
...@@ -497,7 +468,7 @@ void testVdwTaper( FILE* log ) { ...@@ -497,7 +468,7 @@ void testVdwTaper( FILE* log ) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy, log ); setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy );
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 3.5478444e+00; double expectedEnergy = 3.5478444e+00;
...@@ -512,12 +483,12 @@ void testVdwTaper( FILE* log ) { ...@@ -512,12 +483,12 @@ void testVdwTaper( FILE* log ) {
expectedForces[7] = Vec3( -0.0000000e+00, -0.0000000e+00, -0.0000000e+00 ); expectedForces[7] = Vec3( -0.0000000e+00, -0.0000000e+00, -0.0000000e+00 );
double tolerance = 1.0e-04; double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
// test PBC // test PBC
void testVdwPBC( FILE* log ) { void testVdwPBC( ) {
std::string testName = "testVdwPBC"; std::string testName = "testVdwPBC";
...@@ -527,7 +498,7 @@ void testVdwPBC( FILE* log ) { ...@@ -527,7 +498,7 @@ void testVdwPBC( FILE* log ) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy, log ); setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy );
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 1.4949141e+01; double expectedEnergy = 1.4949141e+01;
...@@ -545,14 +516,14 @@ void testVdwPBC( FILE* log ) { ...@@ -545,14 +516,14 @@ void testVdwPBC( FILE* log ) {
// if tapering turned off, then absolute difference < 2.0e-05 // if tapering turned off, then absolute difference < 2.0e-05
double tolerance = 5.0e-04; double tolerance = 5.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
} }
// create box of 216 water molecules // create box of 216 water molecules
void setupAndGetForcesEnergyVdwWater( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff, void setupAndGetForcesEnergyVdwWater( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff,
double boxDimension, int includeVdwDispersionCorrection, double boxDimension, int includeVdwDispersionCorrection,
std::vector<Vec3>& forces, double& energy, FILE* log ){ std::vector<Vec3>& forces, double& energy ){
// beginning of Vdw setup // beginning of Vdw setup
...@@ -1271,7 +1242,7 @@ void setupAndGetForcesEnergyVdwWater( const std::string& sigmaCombiningRule, con ...@@ -1271,7 +1242,7 @@ void setupAndGetForcesEnergyVdwWater( const std::string& sigmaCombiningRule, con
// test employing box of 216 water molecules w/ and w/o dispersion correction // test employing box of 216 water molecules w/ and w/o dispersion correction
void testVdwWater( int includeVdwDispersionCorrection, FILE* log ) { void testVdwWater( int includeVdwDispersionCorrection ) {
std::string testName; std::string testName;
...@@ -1287,7 +1258,7 @@ void testVdwWater( int includeVdwDispersionCorrection, FILE* log ) { ...@@ -1287,7 +1258,7 @@ void testVdwWater( int includeVdwDispersionCorrection, FILE* log ) {
std::vector<Vec3> forces; std::vector<Vec3> forces;
double energy; double energy;
setupAndGetForcesEnergyVdwWater( "CUBIC-MEAN", "HHG", cutoff, boxDimension, includeVdwDispersionCorrection, forces, energy, log ); setupAndGetForcesEnergyVdwWater( "CUBIC-MEAN", "HHG", cutoff, boxDimension, includeVdwDispersionCorrection, forces, energy );
std::vector<Vec3> expectedForces(numberOfParticles); std::vector<Vec3> expectedForces(numberOfParticles);
// initialize expected energy and forces // initialize expected energy and forces
...@@ -1952,7 +1923,7 @@ void testVdwWater( int includeVdwDispersionCorrection, FILE* log ) { ...@@ -1952,7 +1923,7 @@ void testVdwWater( int includeVdwDispersionCorrection, FILE* log ) {
// if tapering turned off, then absolute difference < 2.0e-05 // if tapering turned off, then absolute difference < 2.0e-05
double tolerance = 5.0e-04; double tolerance = 5.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log ); compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance );
// test sigma/epsilon rules for dispersion correction // test sigma/epsilon rules for dispersion correction
...@@ -1975,7 +1946,7 @@ void testVdwWater( int includeVdwDispersionCorrection, FILE* log ) { ...@@ -1975,7 +1946,7 @@ void testVdwWater( int includeVdwDispersionCorrection, FILE* log ) {
expectedEnergies.push_back( 3.2774624e+03 ); expectedEnergies.push_back( 3.2774624e+03 );
for( unsigned int ii = 0; ii < sigmaRules.size(); ii++ ){ for( unsigned int ii = 0; ii < sigmaRules.size(); ii++ ){
setupAndGetForcesEnergyVdwWater( sigmaRules[ii], epsilonRules[ii], cutoff, boxDimension, includeVdwDispersionCorrection, forces, energy, log ); setupAndGetForcesEnergyVdwWater( sigmaRules[ii], epsilonRules[ii], cutoff, boxDimension, includeVdwDispersionCorrection, forces, energy );
testName = "testVdwWaterWithDispersionCorrection_" + sigmaRules[ii] + '_' + epsilonRules[ii]; testName = "testVdwWaterWithDispersionCorrection_" + sigmaRules[ii] + '_' + epsilonRules[ii];
ASSERT_EQUAL_TOL_MOD( expectedEnergies[ii], energy, tolerance, testName ); ASSERT_EQUAL_TOL_MOD( expectedEnergies[ii], energy, tolerance, testName );
} }
...@@ -2045,48 +2016,47 @@ int main(int argc, char* argv[]) { ...@@ -2045,48 +2016,47 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL;
testVdw( log ); testVdw();
// tests using two ammonia molecules // tests using two ammonia molecules
// test VDW w/ sigmaRule=CubicMean and epsilonRule=HHG // test VDW w/ sigmaRule=CubicMean and epsilonRule=HHG
testVdwAmmoniaCubicMeanHhg( log ); testVdwAmmoniaCubicMeanHhg();
// test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic // test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
testVdwAmmoniaArithmeticArithmetic( log ); testVdwAmmoniaArithmeticArithmetic();
// test VDW w/ sigmaRule=Geometric and epsilonRule=Geometric // test VDW w/ sigmaRule=Geometric and epsilonRule=Geometric
testVdwAmmoniaGeometricGeometric( log ); testVdwAmmoniaGeometricGeometric();
// test VDW w/ sigmaRule=CubicMean and epsilonRule=Harmonic // test VDW w/ sigmaRule=CubicMean and epsilonRule=Harmonic
testVdwAmmoniaCubicMeanHarmonic( log ); testVdwAmmoniaCubicMeanHarmonic();
// test w/ cutoff=0.25 nm; single ixn between two particles (0 and 6); force nonzero on // test w/ cutoff=0.25 nm; single ixn between two particles (0 and 6); force nonzero on
// particle 4 due to reduction applied to NH // particle 4 due to reduction applied to NH
// the distance between 0 and 6 is ~ 0.235 so the ixn is in the tapered region // the distance between 0 and 6 is ~ 0.235 so the ixn is in the tapered region
testVdwTaper( log ); testVdwTaper();
// test PBC // test PBC
testVdwPBC( log ); testVdwPBC();
// tests based on box of water // tests based on box of water
int includeVdwDispersionCorrection = 0; int includeVdwDispersionCorrection = 0;
testVdwWater( includeVdwDispersionCorrection, log ); testVdwWater( includeVdwDispersionCorrection);
// includes tests for various combinations of sigma/epsilon rules // includes tests for various combinations of sigma/epsilon rules
// when computing vdw dispersion correction // when computing vdw dispersion correction
includeVdwDispersionCorrection = 1; includeVdwDispersionCorrection = 1;
testVdwWater( includeVdwDispersionCorrection, log ); testVdwWater( includeVdwDispersionCorrection);
// test triclinic boxes // test triclinic boxes
......
plugins/amoeba/platforms/cuda/tests/TestCudaWcaDispersionForce.cpp
View file @ a568bb12
...@@ -56,19 +56,7 @@ extern "C" void registerAmoebaCudaKernelFactories(); ...@@ -56,19 +56,7 @@ extern "C" void registerAmoebaCudaKernelFactories();
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) {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e\n", testName.c_str(), expectedEnergy, energy );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ 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 );
} }
...@@ -77,7 +65,7 @@ void compareForcesEnergy( std::string& testName, double expectedEnergy, double e ...@@ -77,7 +65,7 @@ void compareForcesEnergy( std::string& testName, double expectedEnergy, double e
// test Wca dispersion // test Wca dispersion
void testWcaDispersionAmmonia( FILE* log ) { void testWcaDispersionAmmonia() {
std::string testName = "testWcaDispersionAmmonia"; std::string testName = "testWcaDispersionAmmonia";
...@@ -151,7 +139,7 @@ void testWcaDispersionAmmonia( FILE* log ) { ...@@ -151,7 +139,7 @@ void testWcaDispersionAmmonia( FILE* log ) {
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 );
// Try changing the particle parameters and make sure it's still correct. // Try changing the particle parameters and make sure it's still correct.
...@@ -168,7 +156,7 @@ void testWcaDispersionAmmonia( FILE* log ) { ...@@ -168,7 +156,7 @@ void testWcaDispersionAmmonia( FILE* log ) {
bool exceptionThrown = false; bool exceptionThrown = false;
try { try {
// This should throw an exception. // This should throw an exception.
compareForcesEnergy(testName, state1.getPotentialEnergy(), state2.getPotentialEnergy(), state1.getForces(), state2.getForces(), tolerance, log); compareForcesEnergy(testName, state1.getPotentialEnergy(), state2.getPotentialEnergy(), state1.getForces(), state2.getForces(), tolerance);
} }
catch (std::exception ex) { catch (std::exception ex) {
exceptionThrown = true; exceptionThrown = true;
...@@ -176,7 +164,7 @@ void testWcaDispersionAmmonia( FILE* log ) { ...@@ -176,7 +164,7 @@ void testWcaDispersionAmmonia( FILE* log ) {
ASSERT(exceptionThrown); ASSERT(exceptionThrown);
amoebaWcaDispersionForce->updateParametersInContext(context); amoebaWcaDispersionForce->updateParametersInContext(context);
state1 = context.getState(State::Forces | State::Energy); state1 = context.getState(State::Forces | State::Energy);
compareForcesEnergy(testName, state1.getPotentialEnergy(), state2.getPotentialEnergy(), state1.getForces(), state2.getForces(), tolerance, log); compareForcesEnergy(testName, state1.getPotentialEnergy(), state2.getPotentialEnergy(), state1.getForces(), state2.getForces(), tolerance);
} }
int main(int argc, char* argv[]) { int main(int argc, char* argv[]) {
...@@ -185,11 +173,10 @@ int main(int argc, char* argv[]) { ...@@ -185,11 +173,10 @@ int main(int argc, char* argv[]) {
registerAmoebaCudaKernelFactories(); registerAmoebaCudaKernelFactories();
if (argc > 1) if (argc > 1)
Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1])); Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
FILE* log = NULL;
// test Wca dispersion force using two ammonia molecules // test Wca dispersion force using two ammonia molecules
testWcaDispersionAmmonia( log ); testWcaDispersionAmmonia();
} catch(const std::exception& e) { } catch(const std::exception& e) {
......
plugins/amoeba/serialization/tests/TestSerializeAmoebaGeneralizedKirkwoodForce.cpp
View file @ a568bb12
...@@ -60,14 +60,6 @@ void testSerialization() { ...@@ -60,14 +60,6 @@ void testSerialization() {
stringstream buffer; stringstream buffer;
XmlSerializer::serialize<AmoebaGeneralizedKirkwoodForce>(&force1, "Force", buffer); XmlSerializer::serialize<AmoebaGeneralizedKirkwoodForce>(&force1, "Force", buffer);
#ifdef AMOEBA_DEBUG
if( 0 ){
FILE* filePtr = fopen("GeneralizedKirkwood.xml", "w" );
(void) fprintf( filePtr, "%s", buffer.str().c_str() );
(void) fclose( filePtr );
}
#endif
AmoebaGeneralizedKirkwoodForce* copy = XmlSerializer::deserialize<AmoebaGeneralizedKirkwoodForce>(buffer); AmoebaGeneralizedKirkwoodForce* copy = XmlSerializer::deserialize<AmoebaGeneralizedKirkwoodForce>(buffer);
// Compare the two forces to see if they are identical. // Compare the two forces to see if they are identical.
......
plugins/amoeba/serialization/tests/TestSerializeAmoebaMultipoleForce.cpp
View file @ a568bb12
...@@ -99,13 +99,6 @@ void testSerialization() { ...@@ -99,13 +99,6 @@ void testSerialization() {
stringstream buffer; stringstream buffer;
XmlSerializer::serialize<AmoebaMultipoleForce>(&force1, "Force", buffer); XmlSerializer::serialize<AmoebaMultipoleForce>(&force1, "Force", buffer);
#ifdef AMOEBA_DEBUG
if( 0 ){
FILE* filePtr = fopen("Multipole.xml", "w" );
(void) fprintf( filePtr, "%s", buffer.str().c_str() );
(void) fclose( filePtr );
}
#endif
AmoebaMultipoleForce* copy = XmlSerializer::deserialize<AmoebaMultipoleForce>(buffer); AmoebaMultipoleForce* copy = XmlSerializer::deserialize<AmoebaMultipoleForce>(buffer);
......
plugins/amoeba/serialization/tests/TestSerializeAmoebaTorsionTorsionForce.cpp
View file @ a568bb12
...@@ -94,14 +94,6 @@ void testSerialization() { ...@@ -94,14 +94,6 @@ void testSerialization() {
stringstream buffer; stringstream buffer;
XmlSerializer::serialize<AmoebaTorsionTorsionForce>(&force1, "Force", buffer); XmlSerializer::serialize<AmoebaTorsionTorsionForce>(&force1, "Force", buffer);
#ifdef AMOEBA_DEBUG
if( 0 ){
FILE* filePtr = fopen("TorsionTorsion.xml", "w" );
(void) fprintf( filePtr, "%s", buffer.str().c_str() );
(void) fclose( filePtr );
}
#endif
AmoebaTorsionTorsionForce* copy = XmlSerializer::deserialize<AmoebaTorsionTorsionForce>(buffer); AmoebaTorsionTorsionForce* copy = XmlSerializer::deserialize<AmoebaTorsionTorsionForce>(buffer);
// Compare the two force1s to see if they are identical. // Compare the two force1s to see if they are identical.
......
plugins/amoeba/serialization/tests/TestSerializeAmoebaVdwForce.cpp
View file @ a568bb12
...@@ -65,14 +65,6 @@ void testSerialization() { ...@@ -65,14 +65,6 @@ void testSerialization() {
stringstream buffer; stringstream buffer;
XmlSerializer::serialize<AmoebaVdwForce>(&force1, "Force", buffer); XmlSerializer::serialize<AmoebaVdwForce>(&force1, "Force", buffer);
#ifdef AMOEBA_DEBUG
if( 0 ){
FILE* filePtr = fopen("Vdw.xml", "w" );
(void) fprintf( filePtr, "%s", buffer.str().c_str() );
(void) fclose( filePtr );
}
#endif
AmoebaVdwForce* copy = XmlSerializer::deserialize<AmoebaVdwForce>(buffer); AmoebaVdwForce* copy = XmlSerializer::deserialize<AmoebaVdwForce>(buffer);
// Compare the two forces to see if they are identical. // Compare the two forces to see if they are identical.
......
plugins/amoeba/serialization/tests/TestSerializeAmoebaWcaDispersionForce.cpp
View file @ a568bb12
...@@ -62,14 +62,6 @@ void testSerialization() { ...@@ -62,14 +62,6 @@ void testSerialization() {
stringstream buffer; stringstream buffer;
XmlSerializer::serialize<AmoebaWcaDispersionForce>(&force1, "Force", buffer); XmlSerializer::serialize<AmoebaWcaDispersionForce>(&force1, "Force", buffer);
#ifdef AMOEBA_DEBUG
if( 0 ){
FILE* filePtr = fopen("WcaDispersion.xml", "w" );
(void) fprintf( filePtr, "%s", buffer.str().c_str() );
(void) fclose( filePtr );
}
#endif
AmoebaWcaDispersionForce* copy = XmlSerializer::deserialize<AmoebaWcaDispersionForce>(buffer); AmoebaWcaDispersionForce* copy = XmlSerializer::deserialize<AmoebaWcaDispersionForce>(buffer);
// Compare the two forces to see if they are identical. // Compare the two forces to see if they are identical.
......
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