Commit ce6b5496 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Mods for gcc 4.4.x

parent 626a504d
...@@ -42,7 +42,7 @@ FOREACH(TEST_PROG ${TEST_PROGS}) ...@@ -42,7 +42,7 @@ FOREACH(TEST_PROG ${TEST_PROGS})
ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG}) ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
#TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_CUDA_TARGET} ${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET}) #TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_CUDA_TARGET} ${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET})
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_CUDA_TARGET}) TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_OPENMM__AMOEBA_TARGET} )
ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT}) ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
# Link with static library # Link with static library
......
...@@ -85,11 +85,12 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou ...@@ -85,11 +85,12 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou
} else { } else {
angle = RADIAN*acos(cosine); angle = RADIAN*acos(cosine);
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "getPrefactorsGivenAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealAngle ); (void) fprintf( log, "getPrefactorsGivenAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealAngle );
(void) fflush( log ); (void) fflush( log );
} }
#endif
double deltaIdeal = angle - idealAngle; double deltaIdeal = angle - idealAngle;
double deltaIdeal2 = deltaIdeal*deltaIdeal; double deltaIdeal2 = deltaIdeal*deltaIdeal;
...@@ -128,12 +129,13 @@ static void computeAmoebaHarmonicAngleForce(int bondIndex, std::vector<Vec3>& p ...@@ -128,12 +129,13 @@ static void computeAmoebaHarmonicAngleForce(int bondIndex, std::vector<Vec3>& p
double quarticK = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAngleQuartic(); double quarticK = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAngleQuartic();
double penticK = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAnglePentic(); double penticK = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAnglePentic();
double sexticK = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAngleSextic(); double sexticK = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAngleSextic();
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaHarmonicAngleForce: bond %d [%d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n", (void) fprintf( log, "computeAmoebaHarmonicAngleForce: 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 ); bondIndex, particle1, particle2, particle3, idealAngle, quadraticK, cubicK, quarticK, penticK, sexticK );
(void) fflush( log ); (void) fflush( log );
} }
#endif
double deltaR[2][3]; double deltaR[2][3];
double r2_0 = 0.0; double r2_0 = 0.0;
...@@ -157,10 +159,12 @@ static void computeAmoebaHarmonicAngleForce(int bondIndex, std::vector<Vec3>& p ...@@ -157,10 +159,12 @@ static void computeAmoebaHarmonicAngleForce(int bondIndex, std::vector<Vec3>& p
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 ){ if( log ){
(void) fprintf( log, "dot=%10.3e r2_0=%10.3e r2_1=%10.3e\n", dot, r2_0, r2_1 ); (void) fprintf( log, "dot=%10.3e r2_0=%10.3e r2_1=%10.3e\n", dot, r2_0, r2_1 );
(void) fflush( log ); (void) fflush( log );
} }
#endif
double dEdR; double dEdR;
double energyTerm; double energyTerm;
...@@ -214,6 +218,7 @@ static void computeAmoebaHarmonicAngleForces( Context& context, AmoebaHarmonicAn ...@@ -214,6 +218,7 @@ static void computeAmoebaHarmonicAngleForces( Context& context, AmoebaHarmonicAn
computeAmoebaHarmonicAngleForce(ii, positions, amoebaHarmonicAngleForce, expectedForces, expectedEnergy, log ); computeAmoebaHarmonicAngleForce(ii, positions, amoebaHarmonicAngleForce, expectedForces, expectedEnergy, log );
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaHarmonicAngleForces: expected energy=%14.7e\n", *expectedEnergy ); (void) fprintf( log, "computeAmoebaHarmonicAngleForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){ for( unsigned int ii = 0; ii < positions.size(); ii++ ){
...@@ -221,6 +226,7 @@ static void computeAmoebaHarmonicAngleForces( Context& context, AmoebaHarmonicAn ...@@ -221,6 +226,7 @@ static void computeAmoebaHarmonicAngleForces( Context& context, AmoebaHarmonicAn
} }
(void) fflush( log ); (void) fflush( log );
} }
#endif
return; return;
} }
...@@ -235,6 +241,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicAngleFor ...@@ -235,6 +241,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicAngleFor
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 ){ if( log ){
(void) fprintf( log, "computeAmoebaHarmonicAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() ); (void) fprintf( log, "computeAmoebaHarmonicAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
...@@ -243,6 +250,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicAngleFor ...@@ -243,6 +250,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicAngleFor
} }
(void) fflush( log ); (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 );
...@@ -299,8 +307,10 @@ int main( int numberOfArguments, char* argv[] ) { ...@@ -299,8 +307,10 @@ int main( int numberOfArguments, char* argv[] ) {
//FILE* log = stderr; //FILE* log = stderr;
testOneAngle( log ); testOneAngle( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr ) if( log && log != stderr )
(void) fclose( log ); (void) fclose( log );
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
......
...@@ -100,7 +100,7 @@ static void computeAmoebaHarmonicBondForces( Context& context, AmoebaHarmonicBon ...@@ -100,7 +100,7 @@ static void computeAmoebaHarmonicBondForces( Context& context, AmoebaHarmonicBon
for( int ii = 0; ii < amoebaHarmonicBondForce.getNumBonds(); ii++ ){ for( int ii = 0; ii < amoebaHarmonicBondForce.getNumBonds(); ii++ ){
computeAmoebaHarmonicBondForce(ii, positions, amoebaHarmonicBondForce, expectedForces, expectedEnergy ); computeAmoebaHarmonicBondForce(ii, positions, amoebaHarmonicBondForce, expectedForces, expectedEnergy );
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaHarmonicBondForces: expected energy=%15.7e\n", *expectedEnergy ); (void) fprintf( log, "computeAmoebaHarmonicBondForces: expected energy=%15.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){ for( unsigned int ii = 0; ii < positions.size(); ii++ ){
...@@ -108,6 +108,7 @@ static void computeAmoebaHarmonicBondForces( Context& context, AmoebaHarmonicBon ...@@ -108,6 +108,7 @@ static void computeAmoebaHarmonicBondForces( Context& context, AmoebaHarmonicBon
} }
(void) fflush( log ); (void) fflush( log );
} }
#endif
return; return;
} }
...@@ -120,7 +121,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicBondForc ...@@ -120,7 +121,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicBondForc
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 ){ if( log ){
(void) fprintf( log, "computeAmoebaHarmonicBondForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy() ); (void) fprintf( log, "computeAmoebaHarmonicBondForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
...@@ -129,6 +130,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicBondForc ...@@ -129,6 +130,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicBondForc
} }
(void) fflush( log ); (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 );
...@@ -209,9 +211,10 @@ int main( int numberOfArguments, char* argv[] ) { ...@@ -209,9 +211,10 @@ int main( int numberOfArguments, char* argv[] ) {
//testOneBond( log ); //testOneBond( log );
testTwoBond( log ); testTwoBond( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr ) if( log && log != stderr )
(void) fclose( log ); (void) fclose( log );
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
......
...@@ -87,11 +87,12 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP ...@@ -87,11 +87,12 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP
} else { } else {
angle = RADIAN*acos(cosine); angle = RADIAN*acos(cosine);
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "getPrefactorsGivenInPlaneAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealInPlaneAngle ); (void) fprintf( log, "getPrefactorsGivenInPlaneAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealInPlaneAngle );
(void) fflush( log ); (void) fflush( log );
} }
#endif
double deltaIdeal = angle - idealInPlaneAngle; double deltaIdeal = angle - idealInPlaneAngle;
double deltaIdeal2 = deltaIdeal*deltaIdeal; double deltaIdeal2 = deltaIdeal*deltaIdeal;
...@@ -130,12 +131,13 @@ static void computeAmoebaHarmonicInPlaneAngleForce(int bondIndex, std::vector<V ...@@ -130,12 +131,13 @@ static void computeAmoebaHarmonicInPlaneAngleForce(int bondIndex, std::vector<V
double quarticK = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAngleQuartic(); double quarticK = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAngleQuartic();
double penticK = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAnglePentic(); double penticK = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAnglePentic();
double sexticK = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAngleSextic(); double sexticK = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAngleSextic();
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForce: bond %d [%d %d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n", (void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForce: 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 ); bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK, cubicK, quarticK, penticK, sexticK );
(void) fflush( log ); (void) fflush( log );
} }
#endif
// T = AD x CD // T = AD x CD
// P = B + T*delta // P = B + T*delta
...@@ -176,10 +178,12 @@ static void computeAmoebaHarmonicInPlaneAngleForce(int bondIndex, std::vector<V ...@@ -176,10 +178,12 @@ static void computeAmoebaHarmonicInPlaneAngleForce(int bondIndex, std::vector<V
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 ){ if( log ){
(void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForce: rAp2 or rCp2 <= 0.0\n" ); (void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForce: rAp2 or rCp2 <= 0.0\n" );
(void) fflush( log ); (void) fflush( log );
} }
#endif
return; return;
} }
...@@ -290,7 +294,7 @@ static void computeAmoebaHarmonicInPlaneAngleForces( Context& context, AmoebaHar ...@@ -290,7 +294,7 @@ static void computeAmoebaHarmonicInPlaneAngleForces( Context& context, AmoebaHar
for( int ii = 0; ii < amoebaHarmonicInPlaneAngleForce.getNumAngles(); ii++ ){ for( int ii = 0; ii < amoebaHarmonicInPlaneAngleForce.getNumAngles(); ii++ ){
computeAmoebaHarmonicInPlaneAngleForce(ii, positions, amoebaHarmonicInPlaneAngleForce, expectedForces, expectedEnergy, log ); computeAmoebaHarmonicInPlaneAngleForce(ii, positions, amoebaHarmonicInPlaneAngleForce, expectedForces, expectedEnergy, log );
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForces: expected energy=%14.7e\n", *expectedEnergy ); (void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){ for( unsigned int ii = 0; ii < positions.size(); ii++ ){
...@@ -298,6 +302,7 @@ static void computeAmoebaHarmonicInPlaneAngleForces( Context& context, AmoebaHar ...@@ -298,6 +302,7 @@ static void computeAmoebaHarmonicInPlaneAngleForces( Context& context, AmoebaHar
} }
(void) fflush( log ); (void) fflush( log );
} }
#endif
return; return;
} }
...@@ -311,7 +316,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicInPlaneA ...@@ -311,7 +316,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicInPlaneA
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 ){ if( log ){
(void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() ); (void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
...@@ -320,6 +325,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicInPlaneA ...@@ -320,6 +325,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicInPlaneA
} }
(void) fflush( log ); (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 );
...@@ -377,8 +383,10 @@ int main( int numberOfArguments, char* argv[] ) { ...@@ -377,8 +383,10 @@ int main( int numberOfArguments, char* argv[] ) {
//FILE* log = fopen( "AmoebaHarmonicInPlaneAngleForce.log", "w" );; //FILE* log = fopen( "AmoebaHarmonicInPlaneAngleForce.log", "w" );;
testOneAngle( NULL ); testOneAngle( NULL );
#ifdef AMOEBA_DEBUG
if( log && log != stderr ) if( log && log != stderr )
(void) fclose( log ); (void) fclose( log );
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
......
...@@ -88,12 +88,14 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -88,12 +88,14 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
int particle1, particle2, particle3, particle4; int particle1, particle2, particle3, particle4;
double kAngleQuadratic; double kAngleQuadratic;
amoebaOutOfPlaneBendForce.getOutOfPlaneBendParameters(bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic ); amoebaOutOfPlaneBendForce.getOutOfPlaneBendParameters(bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic );
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForce: bond %d [%d %d %d %d] k=[%10.3e %10.3e %10.3e %10.3e %10.3e]\n", (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 ); bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic, kAngleCubic, kAngleQuartic, kAnglePentic, kAngleSextic );
(void) fflush( log ); (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 };
...@@ -135,12 +137,13 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& ...@@ -135,12 +137,13 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
} else { } else {
angle = RADIAN*acos(cosine); angle = RADIAN*acos(cosine);
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForce: bkk2=%14.7e rDB2=%14.7e cos=%14.7e dt=%14.7e]\n", (void) fprintf( log, "computeAmoebaOutOfPlaneBendForce: bkk2=%14.7e rDB2=%14.7e cos=%14.7e dt=%14.7e]\n",
bkk2, rDB2, cosine, angle ); bkk2, rDB2, cosine, angle );
(void) fflush( log ); (void) fflush( log );
} }
#endif
// chain rule // chain rule
...@@ -256,6 +259,7 @@ static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlan ...@@ -256,6 +259,7 @@ static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlan
computeAmoebaOutOfPlaneBendForce(ii, positions, amoebaOutOfPlaneBendForce, expectedForces, expectedEnergy, log ); computeAmoebaOutOfPlaneBendForce(ii, positions, amoebaOutOfPlaneBendForce, expectedForces, expectedEnergy, log );
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e\n", *expectedEnergy ); (void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){ for( unsigned int ii = 0; ii < positions.size(); ii++ ){
...@@ -263,6 +267,7 @@ static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlan ...@@ -263,6 +267,7 @@ static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlan
} }
(void) fflush( log ); (void) fflush( log );
} }
#endif
return; return;
} }
...@@ -276,7 +281,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaOutOfPlaneBendFo ...@@ -276,7 +281,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaOutOfPlaneBendFo
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 ){ if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() ); (void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
...@@ -285,6 +290,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaOutOfPlaneBendFo ...@@ -285,6 +290,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaOutOfPlaneBendFo
} }
(void) fflush( log ); (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 );
...@@ -412,16 +418,20 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -412,16 +418,20 @@ 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 ){ if( log ){
(void) fprintf( log, "Set id %d not recognized.\n", setId ); (void) fprintf( log, "Set id %d not recognized.\n", setId );
} }
char buffer[1024]; #endif
(void) sprintf( buffer, "Set id %d not recognized.", setId ); std::stringstream buffer;
throw OpenMMException( buffer ); buffer << "Set id " << setId << " not recognized.";
throw OpenMMException( buffer.str() );
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "Set id %d.\n", setId ); (void) fprintf( log, "Set id %d.\n", setId );
} }
#endif
amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend ); amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend );
...@@ -431,12 +441,14 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -431,12 +441,14 @@ 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 ){ if( log ){
(void) fprintf( log, "Coordinates %d not loaded.", particleIndices[ii] ); (void) fprintf( log, "Coordinates %d not loaded.", particleIndices[ii] );
} }
char buffer[1024]; #endif
(void) sprintf( buffer, "Coordinates %d not loaded.", particleIndices[ii] ); std::stringstream buffer;
throw OpenMMException( buffer ); buffer << "Coordinates " << particleIndices[ii] << " not loaded.";
throw OpenMMException( buffer.str() );
} }
positions[ii] = coordinates[particleIndices[ii]]; positions[ii] = coordinates[particleIndices[ii]];
} }
...@@ -475,6 +487,7 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -475,6 +487,7 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
} }
if( iter == 6 ){ if( iter == 6 ){
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: energy=%14.7e\n", totalEnergy); (void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: energy=%14.7e\n", totalEnergy);
for( std::map<int,Vec3>::iterator ii = totalForces.begin(); ii != totalForces.end(); ii++ ){ for( std::map<int,Vec3>::iterator ii = totalForces.begin(); ii != totalForces.end(); ii++ ){
...@@ -485,6 +498,7 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) { ...@@ -485,6 +498,7 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
} }
(void) fflush( log ); (void) fflush( log );
} }
#endif
} }
} }
...@@ -504,8 +518,10 @@ int main( int numberOfArguments, char* argv[] ) { ...@@ -504,8 +518,10 @@ int main( int numberOfArguments, char* argv[] ) {
//for( int ii = 1; ii <= 6; ii++ ){ //for( int ii = 1; ii <= 6; ii++ ){
// testOneOutOfPlaneBend2( log, ii ); // testOneOutOfPlaneBend2( log, ii );
//} //}
#ifdef AMOEBA_DEBUG
if( log && log != stderr ) if( log && log != stderr )
(void) fclose( log ); (void) fclose( log );
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
......
...@@ -82,12 +82,13 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit ...@@ -82,12 +82,13 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
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 ){ if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForce: bond %d [%d %d %d %d %d %d] k=%10.3e\n", (void) fprintf( log, "computeAmoebaPiTorsionForce: bond %d [%d %d %d %d %d %d] k=%10.3e\n",
bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion ); bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion );
(void) fflush( log ); (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];
...@@ -227,7 +228,7 @@ static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce ...@@ -227,7 +228,7 @@ static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce
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, log );
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e\n", *expectedEnergy ); (void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){ for( unsigned int ii = 0; ii < positions.size(); ii++ ){
...@@ -235,6 +236,7 @@ static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce ...@@ -235,6 +236,7 @@ static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce
} }
(void) fflush( log ); (void) fflush( log );
} }
#endif
return; return;
} }
...@@ -248,7 +250,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce& ...@@ -248,7 +250,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce&
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 ){ if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() ); (void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
...@@ -257,6 +259,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce& ...@@ -257,6 +259,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce&
} }
(void) fflush( log ); (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 );
...@@ -307,8 +310,10 @@ int main( int numberOfArguments, char* argv[] ) { ...@@ -307,8 +310,10 @@ int main( int numberOfArguments, char* argv[] ) {
//FILE* log = fopen( "AmoebaPiTorsionForce1.log", "w" );; //FILE* log = fopen( "AmoebaPiTorsionForce1.log", "w" );;
testOnePiTorsion( log ); testOnePiTorsion( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr ) if( log && log != stderr )
(void) fclose( log ); (void) fclose( log );
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
......
...@@ -85,11 +85,13 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -85,11 +85,13 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend); amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend);
angleStretchBend *= RADIAN; angleStretchBend *= RADIAN;
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForce: bond %d [%d %d %d] ab=%10.3e cb=%10.3e angle=%10.3e k=%10.3e\n", (void) fprintf( log, "computeAmoebaStretchBendForce: bond %d [%d %d %d] ab=%10.3e cb=%10.3e angle=%10.3e k=%10.3e\n",
bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend ); bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend );
(void) fflush( log ); (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 };
...@@ -182,11 +184,12 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos ...@@ -182,11 +184,12 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
forces[particle3][2] -= subForce[2][2]; forces[particle3][2] -= subForce[2][2];
*energy += term*dt*dr; *energy += term*dt*dr;
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForce: angle=%10.3e dt=%10.3e dr=%10.3e\n", angle, dt, dr ); (void) fprintf( log, "computeAmoebaStretchBendForce: angle=%10.3e dt=%10.3e dr=%10.3e\n", angle, dt, dr );
(void) fflush( log ); (void) fflush( log );
} }
#endif
return; return;
} }
...@@ -210,7 +213,7 @@ static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendF ...@@ -210,7 +213,7 @@ static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendF
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, log );
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForces: expected energy=%14.7e\n", *expectedEnergy ); (void) fprintf( log, "computeAmoebaStretchBendForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){ for( unsigned int ii = 0; ii < positions.size(); ii++ ){
...@@ -218,6 +221,7 @@ static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendF ...@@ -218,6 +221,7 @@ static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendF
} }
(void) fflush( log ); (void) fflush( log );
} }
#endif
return; return;
} }
...@@ -231,7 +235,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaStretchBendForce ...@@ -231,7 +235,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaStretchBendForce
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 ){ if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() ); (void) fprintf( log, "computeAmoebaStretchBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
...@@ -240,6 +244,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaStretchBendForce ...@@ -240,6 +244,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaStretchBendForce
} }
(void) fflush( log ); (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 );
...@@ -291,8 +296,10 @@ int main( int numberOfArguments, char* argv[] ) { ...@@ -291,8 +296,10 @@ int main( int numberOfArguments, char* argv[] ) {
//FILE* log = stderr; //FILE* log = stderr;
//FILE* log = fopen( "AmoebaStretchBendForce1.log", "w" );; //FILE* log = fopen( "AmoebaStretchBendForce1.log", "w" );;
testOneStretchBend( log ); testOneStretchBend( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr ) if( log && log != stderr )
(void) fclose( log ); (void) fclose( log );
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
......
...@@ -95,6 +95,7 @@ static void computeAmoebaTorsionForce(int bondIndex, std::vector<Vec3>& positio ...@@ -95,6 +95,7 @@ static void computeAmoebaTorsionForce(int bondIndex, std::vector<Vec3>& positio
torsions.push_back( torsion2 ); torsions.push_back( torsion2 );
torsions.push_back( torsion3 ); torsions.push_back( torsion3 );
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaTorsionForce: bond %d [%d %d %d %d]\n", (void) fprintf( log, "computeAmoebaTorsionForce: bond %d [%d %d %d %d]\n",
bondIndex, particle1, particle2, particle3, particle4 ); bondIndex, particle1, particle2, particle3, particle4 );
...@@ -103,6 +104,7 @@ static void computeAmoebaTorsionForce(int bondIndex, std::vector<Vec3>& positio ...@@ -103,6 +104,7 @@ static void computeAmoebaTorsionForce(int bondIndex, std::vector<Vec3>& positio
} }
(void) fflush( log ); (void) fflush( log );
} }
#endif
enum { BA, CB, DC, CA, DB, LastDeltaIndex }; enum { BA, CB, DC, CA, DB, LastDeltaIndex };
double deltaR[LastDeltaIndex][3]; double deltaR[LastDeltaIndex][3];
...@@ -251,6 +253,7 @@ static void computeAmoebaTorsionForces( Context& context, AmoebaTorsionForce& am ...@@ -251,6 +253,7 @@ static void computeAmoebaTorsionForces( Context& context, AmoebaTorsionForce& am
computeAmoebaTorsionForce(ii, positions, amoebaTorsionForce, expectedForces, expectedEnergy, log ); computeAmoebaTorsionForce(ii, positions, amoebaTorsionForce, expectedForces, expectedEnergy, log );
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaTorsionForces: expected energy=%14.7e\n", *expectedEnergy ); (void) fprintf( log, "computeAmoebaTorsionForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){ for( unsigned int ii = 0; ii < positions.size(); ii++ ){
...@@ -258,6 +261,8 @@ static void computeAmoebaTorsionForces( Context& context, AmoebaTorsionForce& am ...@@ -258,6 +261,8 @@ static void computeAmoebaTorsionForces( Context& context, AmoebaTorsionForce& am
} }
(void) fflush( log ); (void) fflush( log );
} }
#endif
return; return;
} }
...@@ -272,6 +277,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaTorsionForce& am ...@@ -272,6 +277,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaTorsionForce& am
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 ){ if( log ){
(void) fprintf( log, "computeAmoebaTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() ); (void) fprintf( log, "computeAmoebaTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
...@@ -280,6 +286,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaTorsionForce& am ...@@ -280,6 +286,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaTorsionForce& am
} }
(void) fflush( log ); (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 );
...@@ -336,8 +343,10 @@ int main( int numberOfArguments, char* argv[] ) { ...@@ -336,8 +343,10 @@ int main( int numberOfArguments, char* argv[] ) {
FILE* log = NULL; FILE* log = NULL;
//FILE* log = fopen( "AmoebaTorsionForce.log", "w" );; //FILE* log = fopen( "AmoebaTorsionForce.log", "w" );;
testOneTorsion( log ); testOneTorsion( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr ) if( log && log != stderr )
(void) fclose( log ); (void) fclose( log );
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
......
...@@ -102,6 +102,7 @@ static void computeAmoebaUreyBradleyForces( Context& context, AmoebaUreyBradleyF ...@@ -102,6 +102,7 @@ static void computeAmoebaUreyBradleyForces( Context& context, AmoebaUreyBradleyF
computeAmoebaUreyBradleyForce(ii, positions, amoebaUreyBradleyForce, expectedForces, expectedEnergy ); computeAmoebaUreyBradleyForce(ii, positions, amoebaUreyBradleyForce, expectedForces, expectedEnergy );
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaUreyBradleyForces: expected energy=%15.7e\n", *expectedEnergy ); (void) fprintf( log, "computeAmoebaUreyBradleyForces: expected energy=%15.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){ for( unsigned int ii = 0; ii < positions.size(); ii++ ){
...@@ -109,6 +110,7 @@ static void computeAmoebaUreyBradleyForces( Context& context, AmoebaUreyBradleyF ...@@ -109,6 +110,7 @@ static void computeAmoebaUreyBradleyForces( Context& context, AmoebaUreyBradleyF
} }
(void) fflush( log ); (void) fflush( log );
} }
#endif
return; return;
} }
...@@ -122,6 +124,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaUreyBradleyForce ...@@ -122,6 +124,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaUreyBradleyForce
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 ){ if( log ){
(void) fprintf( log, "computeAmoebaUreyBradleyForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy() ); (void) fprintf( log, "computeAmoebaUreyBradleyForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
...@@ -130,6 +133,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaUreyBradleyForce ...@@ -130,6 +133,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaUreyBradleyForce
} }
(void) fflush( log ); (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 );
...@@ -211,8 +215,10 @@ int main( int numberOfArguments, char* argv[] ) { ...@@ -211,8 +215,10 @@ int main( int numberOfArguments, char* argv[] ) {
//testOneBond( log ); //testOneBond( log );
testTwoBond( log ); testTwoBond( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr ) if( log && log != stderr )
(void) fclose( log ); (void) fclose( log );
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
......
...@@ -133,7 +133,7 @@ void testVdw( FILE* log ) { ...@@ -133,7 +133,7 @@ void testVdw( FILE* log ) {
forces[ii][1] *= conversion; forces[ii][1] *= conversion;
forces[ii][2] *= conversion; forces[ii][2] *= conversion;
} }
#ifdef AMOEBA_DEBUG
if( log ){ if( log ){
(void) fprintf( log, "computeAmoebaVdwForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() ); (void) fprintf( log, "computeAmoebaVdwForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){ for( unsigned int ii = 0; ii < forces.size(); ii++ ){
...@@ -142,6 +142,7 @@ void testVdw( FILE* log ) { ...@@ -142,6 +142,7 @@ void testVdw( FILE* log ) {
} }
(void) fflush( log ); (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++ ){
...@@ -161,8 +162,10 @@ int main( int numberOfArguments, char* argv[] ) { ...@@ -161,8 +162,10 @@ int main( int numberOfArguments, char* argv[] ) {
//FILE* log = stderr; //FILE* log = stderr;
//FILE* log = fopen( "AmoebaVdwForce1.log", "w" );; //FILE* log = fopen( "AmoebaVdwForce1.log", "w" );;
testVdw( log ); testVdw( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr ) if( log && log != stderr )
(void) fclose( log ); (void) fclose( log );
#endif
} }
catch(const std::exception& e) { catch(const std::exception& e) {
......
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