Commit 0304bced authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Fix for gcc 4.4.x

parent ce6b5496
......@@ -24,6 +24,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include <stdio.h>
#include "AmoebaTinkerParameterFile.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/internal/ContextImpl.h"
......
......@@ -87,10 +87,12 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou
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 deltaIdeal2 = deltaIdeal*deltaIdeal;
......@@ -130,11 +132,13 @@ static void computeAmoebaHarmonicAngleForce(int bondIndex, std::vector<Vec3>& p
double penticK = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAnglePentic();
double sexticK = amoebaHarmonicAngleForce.getAmoebaGlobalHarmonicAngleSextic();
#ifdef AMOEBA_DEBUG
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",
bondIndex, particle1, particle2, particle3, idealAngle, quadraticK, cubicK, quarticK, penticK, sexticK );
(void) fflush( log );
}
#endif
double deltaR[2][3];
double r2_0 = 0.0;
......@@ -158,10 +162,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 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 energyTerm;
......@@ -215,6 +221,7 @@ static void computeAmoebaHarmonicAngleForces( Context& context, AmoebaHarmonicAn
computeAmoebaHarmonicAngleForce(ii, positions, amoebaHarmonicAngleForce, expectedForces, expectedEnergy, log );
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaHarmonicAngleForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
......@@ -222,6 +229,8 @@ static void computeAmoebaHarmonicAngleForces( Context& context, AmoebaHarmonicAn
}
(void) fflush( log );
}
#endif
return;
}
......@@ -236,6 +245,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicAngleFor
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaHarmonicAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
......@@ -244,6 +254,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicAngleFor
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
......@@ -299,8 +310,10 @@ int main( int numberOfArguments, char* argv[] ) {
FILE* log = NULL;
testOneAngle( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
#endif
}
catch(const std::exception& e) {
......
......@@ -102,6 +102,7 @@ static void computeAmoebaHarmonicBondForces( Context& context, AmoebaHarmonicBon
computeAmoebaHarmonicBondForce(ii, positions, amoebaHarmonicBondForce, expectedForces, expectedEnergy );
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaHarmonicBondForces: expected energy=%15.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
......@@ -109,6 +110,7 @@ static void computeAmoebaHarmonicBondForces( Context& context, AmoebaHarmonicBon
}
(void) fflush( log );
}
#endif
return;
}
......@@ -122,6 +124,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicBondForc
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaHarmonicBondForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
......@@ -130,6 +133,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicBondForc
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
......@@ -210,9 +214,10 @@ int main( int numberOfArguments, char* argv[] ) {
//testOneBond( log );
testTwoBond( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
#endif
}
catch(const std::exception& e) {
......
......@@ -88,10 +88,12 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP
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 deltaIdeal2 = deltaIdeal*deltaIdeal;
......@@ -131,11 +133,13 @@ static void computeAmoebaHarmonicInPlaneAngleForce(int bondIndex, std::vector<V
double penticK = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAnglePentic();
double sexticK = amoebaHarmonicInPlaneAngleForce.getAmoebaGlobalHarmonicInPlaneAngleSextic();
#ifdef AMOEBA_DEBUG
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",
bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK, cubicK, quarticK, penticK, sexticK );
(void) fflush( log );
}
#endif
// T = AD x CD
// P = B + T*delta
......@@ -176,10 +180,13 @@ static void computeAmoebaHarmonicInPlaneAngleForce(int bondIndex, std::vector<V
double rAp2 = dotVector3( deltaR[AP], deltaR[AP] );
double rCp2 = dotVector3( deltaR[CP], deltaR[CP] );
if( rAp2 <= 0.0 && rCp2 <= 0.0 ){
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForce: rAp2 or rCp2 <= 0.0\n" );
(void) fflush( log );
}
#endif
return;
}
......@@ -291,6 +298,7 @@ static void computeAmoebaHarmonicInPlaneAngleForces( Context& context, AmoebaHar
computeAmoebaHarmonicInPlaneAngleForce(ii, positions, amoebaHarmonicInPlaneAngleForce, expectedForces, expectedEnergy, log );
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
......@@ -298,6 +306,7 @@ static void computeAmoebaHarmonicInPlaneAngleForces( Context& context, AmoebaHar
}
(void) fflush( log );
}
#endif
return;
}
......@@ -312,6 +321,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicInPlaneA
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaHarmonicInPlaneAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
......@@ -320,6 +330,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaHarmonicInPlaneA
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
......@@ -376,8 +387,10 @@ int main( int numberOfArguments, char* argv[] ) {
//FILE* log = fopen( "AmoebaHarmonicInPlaneAngleForce.log", "w" );;
testOneAngle( NULL );
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
#endif
}
catch(const std::exception& e) {
......
......@@ -89,11 +89,13 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
double 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 { AB, CB, DB, AD, CD, LastDeltaIndex };
......@@ -136,11 +138,13 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
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
......@@ -256,6 +260,7 @@ static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlan
computeAmoebaOutOfPlaneBendForce(ii, positions, amoebaOutOfPlaneBendForce, expectedForces, expectedEnergy, log );
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
......@@ -263,6 +268,8 @@ static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlan
}
(void) fflush( log );
}
#endif
return;
}
......@@ -277,6 +284,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaOutOfPlaneBendFo
State state = context.getState(State::Forces | State::Energy);
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++ ){
......@@ -285,6 +293,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaOutOfPlaneBendFo
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
......@@ -412,16 +421,20 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
particleIndices.push_back( 442 );
kOutOfPlaneBend = 0.214755281E-01;
} else {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "Set id %d not recognized.\n", setId );
}
char buffer[1024];
(void) sprintf( buffer, "Set id %d not recognized.", setId );
throw OpenMMException( buffer );
#endif
std::stringstream buffer;
buffer << "Set id " << setId << " not recognized.";
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 );
......@@ -431,12 +444,15 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){
if( coordinates.find( particleIndices[ii] ) == coordinates.end() ){
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "Coordinates %d not loaded.", particleIndices[ii] );
}
char buffer[1024];
(void) sprintf( buffer, "Coordinates %d not loaded.", particleIndices[ii] );
throw OpenMMException( buffer );
#endif
std::stringstream buffer;
buffer << "Coordinates " << particleIndices[ii] << " not loaded.";
throw OpenMMException( buffer.str() );
}
positions[ii] = coordinates[particleIndices[ii]];
}
......@@ -475,6 +491,7 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
}
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++ ){
......@@ -485,6 +502,7 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
}
(void) fflush( log );
}
#endif
}
}
......@@ -504,8 +522,10 @@ int main( int numberOfArguments, char* argv[] ) {
//for( int ii = 1; ii <= 6; ii++ ){
// testOneOutOfPlaneBend2( log, ii );
//}
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
#endif
}
catch(const std::exception& e) {
......
......@@ -109,12 +109,15 @@ void testPMEWater() {
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
#ifdef AMOEBA_DEBUG
(void) fprintf( stderr, "PME forces\n" );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( stderr, "%6u [%14.7e %14.7e %14.7e]\n", ii,
forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( stderr );
#endif
}
......
......@@ -84,11 +84,13 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
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 };
double deltaR[LastDeltaIndex][3];
......@@ -229,6 +231,7 @@ static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce
computeAmoebaPiTorsionForce(ii, positions, amoebaPiTorsionForce, expectedForces, expectedEnergy, log );
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
......@@ -236,6 +239,7 @@ static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce
}
(void) fflush( log );
}
#endif
return;
}
......@@ -250,6 +254,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce&
State state = context.getState(State::Forces | State::Energy);
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++ ){
......@@ -258,6 +263,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce&
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
......@@ -308,8 +314,10 @@ int main( int numberOfArguments, char* argv[] ) {
//FILE* log = fopen( "AmoebaPiTorsionForce1.log", "w" );;
testOnePiTorsion( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
#endif
}
catch(const std::exception& e) {
......
......@@ -84,11 +84,13 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend);
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 k=%10.3e\n",
bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend );
(void) fflush( log );
}
#endif
enum { A, B, C, LastAtomIndex };
enum { AB, CB, CBxAB, ABxP, CBxP, LastDeltaIndex };
......@@ -182,10 +184,12 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
*energy += term*dt*dr;
#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;
}
......@@ -210,6 +214,7 @@ static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendF
computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy, log );
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
......@@ -217,6 +222,7 @@ static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendF
}
(void) fflush( log );
}
#endif
return;
}
......@@ -231,6 +237,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaStretchBendForce
State state = context.getState(State::Forces | State::Energy);
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++ ){
......@@ -239,6 +246,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaStretchBendForce
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
......@@ -290,8 +298,10 @@ int main( int numberOfArguments, char* argv[] ) {
//FILE* log = stderr;
//FILE* log = fopen( "AmoebaStretchBendForce1.log", "w" );;
testOneStretchBend( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
#endif
}
catch(const std::exception& e) {
......
......@@ -95,6 +95,7 @@ static void computeAmoebaTorsionForce(int bondIndex, std::vector<Vec3>& positio
torsions.push_back( torsion2 );
torsions.push_back( torsion3 );
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaTorsionForce: bond %d [%d %d %d %d]\n",
bondIndex, particle1, particle2, particle3, particle4 );
......@@ -103,6 +104,7 @@ static void computeAmoebaTorsionForce(int bondIndex, std::vector<Vec3>& positio
}
(void) fflush( log );
}
#endif
enum { BA, CB, DC, CA, DB, LastDeltaIndex };
double deltaR[LastDeltaIndex][3];
......@@ -251,6 +253,7 @@ static void computeAmoebaTorsionForces( Context& context, AmoebaTorsionForce& am
computeAmoebaTorsionForce(ii, positions, amoebaTorsionForce, expectedForces, expectedEnergy, log );
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaTorsionForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
......@@ -258,6 +261,7 @@ static void computeAmoebaTorsionForces( Context& context, AmoebaTorsionForce& am
}
(void) fflush( log );
}
#endif
return;
}
......@@ -272,6 +276,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaTorsionForce& am
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
......@@ -280,6 +285,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaTorsionForce& am
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
......@@ -336,8 +342,10 @@ int main( int numberOfArguments, char* argv[] ) {
FILE* log = NULL;
//FILE* log = fopen( "AmoebaTorsionForce.log", "w" );;
testOneTorsion( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
#endif
}
catch(const std::exception& e) {
......
......@@ -2660,6 +2660,7 @@ void testTorsionTorsion( FILE* log, int systemId ) {
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++ ){
......@@ -2669,6 +2670,7 @@ void testTorsionTorsion( FILE* log, int systemId ) {
}
(void) fflush( log );
}
#endif
double tolerance = 1.0e-03;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
......@@ -2689,8 +2691,10 @@ int main( int numberOfArguments, char* argv[] ) {
//FILE* log = fopen( "AmoebaTorsionTorsionForce1.log", "w" );;
// testTorsionTorsion( log, 0 );
testTorsionTorsion( log, 1 );
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
#endif
}
catch(const std::exception& e) {
......
......@@ -102,6 +102,7 @@ static void computeAmoebaUreyBradleyForces( Context& context, AmoebaUreyBradleyF
computeAmoebaUreyBradleyForce(ii, positions, amoebaUreyBradleyForce, expectedForces, expectedEnergy );
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaUreyBradleyForces: expected energy=%15.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
......@@ -109,6 +110,7 @@ static void computeAmoebaUreyBradleyForces( Context& context, AmoebaUreyBradleyF
}
(void) fflush( log );
}
#endif
return;
}
......@@ -122,6 +124,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaUreyBradleyForce
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaUreyBradleyForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
......@@ -130,6 +133,7 @@ void compareWithExpectedForceAndEnergy( Context& context, AmoebaUreyBradleyForce
}
(void) fflush( log );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
......@@ -211,8 +215,10 @@ int main( int numberOfArguments, char* argv[] ) {
//testOneBond( log );
testTwoBond( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
#endif
}
catch(const std::exception& e) {
......
......@@ -151,6 +151,7 @@ void testVdw( FILE* log ) {
}
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++ ){
......@@ -159,6 +160,7 @@ void testVdw( FILE* log ) {
}
(void) fflush( log );
}
#endif
double tolerance = 1.0e-03;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
......@@ -178,8 +180,10 @@ int main( int numberOfArguments, char* argv[] ) {
//FILE* log = stderr;
//FILE* log = fopen( "AmoebaVdwForce1.log", "w" );;
testVdw( log );
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
#endif
}
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