Skip to content

GitLab

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

Updated tests w/ energy checks

parent 22ecd657
Hide whitespace changes
Inline Side-by-side
Showing with 482 additions and 63 deletions
platforms/brook/tests/TestBrookStream.cpp
View file @ a01da386
......@@ -47,6 +47,7 @@
#include "OpenMMContext.h"
#include "StandardMMForceField.h"
#include "GBSAOBCForceField.h"
#include "System.h"
#include "LangevinIntegrator.h"
#include "BrookRandomNumberGenerator.h"
......@@ -58,11 +59,13 @@
#include <fstream>
typedef std::vector<std::string> StringVector;
typedef StringVector::iterator StringVectorI;
typedef StringVector::const_iterator StringVectorCI;
using namespace OpenMM;
using namespace std;
const double TOL = 5e-4;
const double TOL = 1e-5;
void testWriteRead( void ){
......@@ -1086,7 +1089,7 @@ void testBrookBonds( void ){
//(void) fprintf( log, "testBrookBonds:Calling getState\n");
//(void) fflush( log );
State state = context.getState( State::Forces );
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
(void) fprintf( log, "Harmonic bond forces\n");
......@@ -1099,9 +1102,10 @@ void testBrookBonds( void ){
ASSERT_EQUAL_VEC(Vec3(0.7*0.2, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
(void) fprintf( log, "Harmonic bond forces ok -- skipping energy test\n");
ASSERT_EQUAL_TOL( 0.5*0.8*0.5*0.5 + 0.5*0.7*0.2*0.2, state.getPotentialEnergy(), TOL);
(void) fprintf( log, "Harmonic bond forces ok\n");
// ASSERT_EQUAL_TOL(0.5*0.8*0.5*0.5 + 0.5*0.7*0.2*0.2, state.getPotentialEnergy(), TOL);
}
void testBrookAngles( void ){
......@@ -1149,7 +1153,7 @@ void testBrookAngles( void ){
//(void) fprintf( log, "%s :Calling getState\n", methodName.c_str() );
//(void) fflush( log );
State state = context.getState( State::Forces );
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
(void) fprintf( log, "Angle bond forces\n");
......@@ -1158,18 +1162,19 @@ void testBrookAngles( void ){
}
(void) fflush( log );
double tolerance = 1.0e-03;
double torque1 = 1.1*PI_M/6;
double torque2 = 1.2*PI_M/4;
ASSERT_EQUAL_VEC(Vec3(torque1, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-0.5*torque2, 0.5*torque2, 0), forces[3], TOL); // reduced by sqrt(2) due to the bond length, another sqrt(2) due to the angle
ASSERT_EQUAL_VEC(Vec3(torque1, 0, 0), forces[0], tolerance);
ASSERT_EQUAL_VEC(Vec3(-0.5*torque2, 0.5*torque2, 0), forces[3], tolerance); // reduced by sqrt(2) due to the bond length, another sqrt(2) due to the angle
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0],
forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1],
forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]),
Vec3(0, 0, 0), TOL);
Vec3(0, 0, 0), tolerance);
(void) fprintf( log, "Angle bond forces ok -- skipping energy test\n");
ASSERT_EQUAL_TOL(0.5*1.1*(PI_M/6)*(PI_M/6) + 0.5*1.2*(PI_M/4)*(PI_M/4), state.getPotentialEnergy(), tolerance);
// ASSERT_EQUAL_TOL(0.5*1.1*(PI_M/6)*(PI_M/6) + 0.5*1.2*(PI_M/4)*(PI_M/4), state.getPotentialEnergy(), TOL);
(void) fprintf( log, "Angle bond forces ok tolerance=%.2e\n", tolerance );
}
void testBrookPeriodicTorsions( void ){
......@@ -1216,7 +1221,7 @@ void testBrookPeriodicTorsions( void ){
//(void) fprintf( log, "%s :Calling getState\n", methodName.c_str() );
//(void) fflush( log );
State state = context.getState( State::Forces );
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
(void) fprintf( log, "Periodic torsion bond forces\n");
......@@ -1226,18 +1231,22 @@ void testBrookPeriodicTorsions( void ){
(void) fflush( log );
double torque = -2*1.1*std::sin(2*PI_M/3);
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], TOL);
double torque = -2*1.1*std::sin(2*PI_M/3);
double tolerance = 1.0e-03;
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], tolerance );
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, 0), forces[3], tolerance );
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0],
forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1],
forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]),
Vec3(0, 0, 0), TOL);
// ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), TOL);
Vec3(0, 0, 0), tolerance );
(void) fprintf( log, "Periodic torsion bond forces ok -- skipping energy test\n");
//double e1 = 1.1*(1+std::cos(2*PI_M/3));
//double e2 = state.getPotentialEnergy();
//(void) fprintf( log, "E1=%7e E2=%.7e\n", e1, e2 );
//(void) fflush( log );
// ASSERT_EQUAL_TOL(0.5*1.1*(PI_M/6)*(PI_M/6) + 0.5*1.2*(PI_M/4)*(PI_M/4), state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), tolerance );
(void) fprintf( log, "Periodic torsion bond forces ok -- tolerance=%.2e\n", tolerance);
}
void testBrookRBTorsions( void ){
......@@ -1264,6 +1273,10 @@ void testBrookRBTorsions( void ){
StandardMMForceField* forceField = new StandardMMForceField( numberOfAtoms, 0, 0, 0, 1 );
for( int ii = 0; ii < numberOfAtoms; ii++ ){
forceField->setAtomParameters( ii, 0.0, 1, 0);
}
forceField->setRBTorsionParameters(0, 0, 1, 2, 3, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6);
system.addForce(forceField);
......@@ -1283,7 +1296,8 @@ void testBrookRBTorsions( void ){
//(void) fprintf( log, "%s :Calling getState\n", methodName.c_str() );
//(void) fflush( log );
State state = context.getState( State::Forces );
//State state = context.getState( State::Forces | State::Energy );
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
(void) fprintf( log, "RB torsion bond forces\n");
......@@ -1298,25 +1312,23 @@ void testBrookRBTorsions( void ){
double c = 0.1*(i+1);
torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
}
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, -0.5*torque), forces[3], TOL);
double tolerance = 0.1;
ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], tolerance );
ASSERT_EQUAL_VEC(Vec3(0, 0.5*torque, -0.5*torque), forces[3], tolerance );
ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0],
forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1],
forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]),
Vec3(0, 0, 0), TOL);
/*
double energy = 0.0;
for (int i = 0; i < 6; ++i) {
double c = 0.1*(i+1);
energy += c*std::pow(std::cos(psi), i);
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
*/
Vec3(0, 0, 0), tolerance);
double energy = 0.0;
for (int i = 0; i < 6; ++i) {
double c = 0.1*(i+1);
energy += c*std::pow(std::cos(psi), i);
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), tolerance );
(void) fprintf( log, "RB torsion bond forces ok -- skipping energy test\n");
// ASSERT_EQUAL_TOL(0.5*1.1*(PI_M/6)*(PI_M/6) + 0.5*1.2*(PI_M/4)*(PI_M/4), state.getPotentialEnergy(), TOL);
(void) fprintf( log, "RB torsion bond forces ok tolerance=%.2e\n", tolerance); fflush( log );
}
void testBrookCoulomb( void ){
......@@ -1360,7 +1372,7 @@ void testBrookCoulomb( void ){
//(void) fprintf( log, "%s :Calling getState\n", methodName.c_str() );
//(void) fflush( log );
State state = context.getState( State::Forces );
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
(void) fprintf( log, "Coulomb forces\n");
......@@ -1372,9 +1384,11 @@ void testBrookCoulomb( void ){
double force = 138.935485*(-0.75)/4.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
//ASSERT_EQUAL_TOL(138.935485*(-0.75)/2.0, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_TOL(138.935485*(-0.75)/2.0, state.getPotentialEnergy(), TOL);
(void) fprintf( log, "Coulomb forces ok -- skipping energy test\n");
(void) fprintf( log, "Coulomb forces ok\n"); fflush( log );
// delete forceField;
}
......@@ -1420,7 +1434,7 @@ void testBrookLJ( void ){
//(void) fprintf( log, "%s :Calling getState\n", methodName.c_str() );
//(void) fflush( log );
State state = context.getState( State::Forces );
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
(void) fprintf( log, "LJ forces\n");
......@@ -1434,9 +1448,9 @@ void testBrookLJ( void ){
double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/2.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
//ASSERT_EQUAL_TOL(4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0)), state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_TOL(4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0)), state.getPotentialEnergy(), TOL);
(void) fprintf( log, "LJ forces ok -- skipping energy test\n");
(void) fprintf( log, "LJ forces ok\n"); fflush( log );
}
......@@ -1492,7 +1506,7 @@ void testBrookExclusionsAnd14( void ){
forceField->setAtomParameters(ii, 0, 1.5, 1);
context.reinitialize();
context.setPositions(positions);
State state = context.getState( State::Forces );
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
double x = 1.5/r;
double eps = 1.0;
......@@ -1515,7 +1529,7 @@ void testBrookExclusionsAnd14( void ){
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[ii], TOL);
(void) fprintf( log, "14 LJ forces ok for index=%d\n\n", ii );
(void) fflush( log );
//ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
// Test Coulomb forces
......@@ -1523,7 +1537,7 @@ void testBrookExclusionsAnd14( void ){
forceField->setAtomParameters(ii, 2, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState( State::Forces );
state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces2 = state.getForces();
force = 138.935485*4/(r*r);
energy = 138.935485*4/r;
......@@ -1544,10 +1558,372 @@ void testBrookExclusionsAnd14( void ){
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[ii], TOL);
(void) fprintf( log, "14 Coulomb forces ok for index=%d\n\n", ii );
(void) fflush( log );
// ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
(void) fprintf( log, "ExclusionsAnd14 ok\n"); fflush( log );
}
static OpenMMContext* testObcForceSetup( int numAtoms, int brookContext ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testObcForceSetup";
static const int debug = 1;
static const int maxPrint = 10;
float epsilon = 1.0e-04f;
FILE* log = stdout;
int totalErrors = 0;
// ---------------------------------------------------------------------------------------
Platform* platform;
if( brookContext ){
platform = new BrookPlatform( 32, "cal", log );
} else {
platform = new ReferencePlatform();
}
System* system = new System( numAtoms, 0 );
LangevinIntegrator* integrator = new LangevinIntegrator(0, 0.1, 0.01);
GBSAOBCForceField* forceField = new GBSAOBCForceField(numAtoms);
for (int i = 0; i < numAtoms; ++i){
// charge radius scalingFactor
forceField->setAtomParameters(i, i%2 == 0 ? -1 : 1, 0.15, 1);
//forceField->setAtomParameters(i, i%2 == 0 ? -1 : 1, 1.5, 1);
}
system->addForce(forceField);
OpenMMContext* context = new OpenMMContext( *system, *integrator, *platform );
// Set random positions for all the atoms.
vector<Vec3> positions(numAtoms);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i){
positions[i] = Vec3(1.0*genrand_real2(), 1.0*genrand_real2(), 1.0*genrand_real2());
}
(void) fprintf( log, "ExclusionsAnd14 ok -- skipping energy test\n");
context->setPositions(positions);
return context;
}
/**---------------------------------------------------------------------------------------
Replacement of sorts for strtok() (static method) (Simbios)
Used to parse parameter file lines
Should be moved to Utilities file
@param lineBuffer string to tokenize
@param delimiter token delimter
@return number of args; if return value equals maxTokens, then more tokens than allocated
--------------------------------------------------------------------------------------- */
static char* localStrsep( char** lineBuffer, const char* delimiter ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "\nTinkerParameterSet::strsep"
char *s;
const char *spanp;
int c, sc;
char *tok;
// ---------------------------------------------------------------------------------------
s = *lineBuffer;
if( s == NULL ){
return (NULL);
}
for( tok = s;; ){
c = *s++;
spanp = delimiter;
do {
if( (sc = *spanp++) == c ){
if( c == 0 ){
s = NULL;
} else {
s[-1] = 0;
}
*lineBuffer = s;
return( tok );
}
} while( sc != 0 );
}
}
static OpenMMContext* testObcForceFileSetup( std::string fileName, int brookContext, int* numAtoms ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testObcForceFileSetup";
static const int debug = 1;
static const int maxPrint = 10;
float epsilon = 1.0e-04f;
FILE* log = stdout;
int totalErrors = 0;
const int bufferSize = 1024;
char buffer[bufferSize];
// ---------------------------------------------------------------------------------------
Platform* platform;
if( brookContext ){
platform = new BrookPlatform( 32, "cal", log );
} else {
platform = new ReferencePlatform();
}
LangevinIntegrator* integrator = new LangevinIntegrator(0, 0.1, 0.01);
FILE* readFile = fopen( fileName.c_str(), "r" );
if( !readFile ){
(void) fprintf( log, "\n%s File=%s not opened.", methodName.c_str(), fileName.c_str() );
(void) fflush( log );
return NULL;
} else {
(void) fprintf( log, "\n%s File=%s opened.", methodName.c_str(), fileName.c_str() );
(void) fflush( log );
}
// atom count
int lineCount = 0;
std::string delimiter = " ";
(void) fgets( buffer, bufferSize, readFile );
StringVector tokens;
tokenizeString( buffer, tokens, delimiter );
if( tokens.size() < 1 ){
std::stringstream message;
message << "\nProblem w/ line=" << lineCount << " <" << buffer << ">";
(void) fprintf( log, "\n%s", message.str().c_str() );
return NULL;
}
StringVectorI ii = tokens.begin();
int numberOfAtoms = atoi( (*ii).c_str() );
(void) fprintf( log, "\n%s %d\n", methodName.c_str(), numberOfAtoms );
*numAtoms = numberOfAtoms;
lineCount++;
System* system = new System( *numAtoms, 0 );
GBSAOBCForceField* forceField = new GBSAOBCForceField(numberOfAtoms);
vector<Vec3> positions(numberOfAtoms);
int index = 0;
for (int i = 0; i < numberOfAtoms; ++i){
(void) fgets( buffer, bufferSize, readFile );
StringVector tokens;
tokenizeString( buffer, tokens, delimiter );
if( tokens.size() < 8 ){
std::stringstream message;
message << "\nProblem w/ line=" << lineCount << " <" << buffer << ">";
(void) fprintf( log, "\n%s", message.str().c_str() );
return NULL;
}
StringVectorI ii = tokens.begin();
ii++;
float coordX = atof( (*ii).c_str() ); ii++;
float coordY = atof( (*ii).c_str() ); ii++;
float coordZ = atof( (*ii).c_str() ); ii++;
float radius = atof( (*ii).c_str() ); ii++;
float scalingFactor = atof( (*ii).c_str() ); ii++;
float charge = atof( (*ii).c_str() ); ii++;
float bornRadi = atof( (*ii).c_str() ); ii++;
positions[index++] = Vec3( coordX, coordY, coordZ );
// charge radius scalingFactor
forceField->setAtomParameters( i, charge, radius, scalingFactor );
(void) fprintf( log, "%d [%.6f %.6f %.6f] q=%.6f rad=%.6f scl=%.6f bR=%.6f\n", i,
coordX, coordY, coordZ, charge, radius, scalingFactor, bornRadi );
lineCount++;
}
(void) fflush( log );
system->addForce(forceField);
OpenMMContext* context = new OpenMMContext( *system, *integrator, *platform );
(void) fflush( log );
// Set positions for all the atoms.
context->setPositions(positions);
return context;
}
void testObcForce() {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testObcForce";
static const int debug = 1;
static const int maxPrint = 10;
float epsilon = 1.0e-04f;
FILE* log = stdout;
int totalErrors = 0;
// ---------------------------------------------------------------------------------------
int numAtoms = 10;
//OpenMMContext* context = testObcForceSetup( numAtoms, 0 );
//OpenMMContext* brookContext = testObcForceSetup( numAtoms, 1 );
OpenMMContext* context = testObcForceFileSetup( std::string( "ObcInfo.txt" ), 0, &numAtoms );
//OpenMMContext* context = NULL;
OpenMMContext* brookContext = testObcForceFileSetup( std::string( "ObcInfo.txt" ), 1, &numAtoms );
//OpenMMContext* brookContext = NULL;
State* state;
vector<Vec3> forces;
if( context ){
State state = context->getState( State::Forces );
forces = state.getForces();
}
vector<Vec3> brookForces;
if( brookContext ){
State brookState = brookContext->getState( State::Forces );
brookForces = brookState.getForces();
}
(void) fprintf( log, "%s OBC forces [%s %s]\n", methodName.c_str(), (context?"Ref":""), (brookContext?"Brook":"") );
for( int ii = 0; ii < numAtoms; ii++ ){
(void) fprintf( log, "%4d ", ii );
if( context && brookContext ){
double diff[3];
double rdiff[3];
int hit = 0;
for( int jj = 0; jj < 3; jj++ ){
diff[jj] = fabs( forces[ii][jj] - brookForces[ii][jj] );
if( forces[ii][jj] != 0.0 ){
rdiff[jj] = fabs( diff[jj]/forces[ii][jj] );
if( rdiff[jj] > 0.01 ){
hit++;
}
} else {
rdiff[jj] = diff[jj];
}
}
(void) fprintf( log, "diff[%8.3f %8.3f %8.3f] [%8.3f %8.3f %8.3f] ", rdiff[0], rdiff[1], rdiff[2], diff[0], diff[1], diff[2] );
if( hit ){
(void) fprintf( log, " XXX" );
}
}
if( context ){
(void) fprintf( log, "[%14.5e %14.5e %14.5e] ", forces[ii][0], forces[ii][1], forces[ii][2] );
}
if( brookContext ){
(void) fprintf( log, "[%14.5e %14.5e %14.5e]", brookForces[ii][0], brookForces[ii][1], brookForces[ii][2] );
}
(void) fprintf( log, "\n" );
}
(void) fflush( log );
double tolerance = 1.0e-03;
if( context && brookContext ){
for (int i = 0; i < numAtoms; ++i) {
Vec3 f = forces[i];
Vec3 fBrook = brookForces[i];
ASSERT_EQUAL_VEC( f, fBrook, tolerance );
}
}
(void) fprintf( log, "testObcForce ok w/ tolerance=%.3e\n", tolerance );
(void) fflush( log );
}
void testObcSingleAtom() {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testObcSingleAtom";
FILE* log = stdout;
// ---------------------------------------------------------------------------------------
BrookPlatform platform;
//ReferencePlatform platform;
System system(1, 0);
system.setAtomMass(0, 2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForceField* forceField = new GBSAOBCForceField(1);
forceField->setAtomParameters(0, 0.5, 0.15, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double bornRadius = 0.15-0.009; // dielectric offset
double eps0 = EPSILON0;
double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric())/bornRadius;
double extendedRadius = bornRadius+0.14; // probe radius
double nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.15/bornRadius, 6.0); // Where did this formula come from? Just copied it from CpuImplicitSolvent.cpp
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
(void) fprintf( log, "%s ok\n", methodName.c_str() );
(void) fflush( log );
}
void testObcEConsistentForce() {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testObcEConsistentForce";
FILE* log = stdout;
// ---------------------------------------------------------------------------------------
BrookPlatform platform;
//ReferencePlatform platform;
const int numAtoms = 10;
System system(numAtoms, 0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForceField* forceField = new GBSAOBCForceField(numAtoms);
for (int i = 0; i < numAtoms; ++i)
forceField->setAtomParameters(i, i%2 == 0 ? -1 : 1, 0.15, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
// Set random positions for all the atoms.
vector<Vec3> positions(numAtoms);
init_gen_rand(0);
for (int i = 0; i < numAtoms; ++i)
positions[i] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
// Take a small step in the direction of the energy gradient.
double norm = 0.0;
for (int i = 0; i < numAtoms; ++i) {
Vec3 f = state.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-3;
double step = delta/norm;
for (int i = 0; i < numAtoms; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
context.setPositions(positions);
// See whether the potential energy changed by the expected amount.
State state2 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta, 0.01)
(void) fprintf( log, "%s ok\n", methodName.c_str() ); (void) fflush( log );
}
int testBrookStreams( int streamSize ){
......@@ -1677,11 +2053,11 @@ int testBrookStreams( int streamSize ){
return totalErrors;
}
void testLangevinSingleBond() {
static OpenMMContext* testLangevinSingleBondSetup( int brookContext, LangevinIntegrator** outIntegrator ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "LangevinSingleBond";
static const std::string methodName = "LangevinSingleBondSetup";
static const int debug = 1;
FILE* log = stdout;
......@@ -1690,41 +2066,77 @@ void testLangevinSingleBond() {
// ---------------------------------------------------------------------------------------
(void) fprintf( log, "%s\n", methodName.c_str() );
(void) fprintf( log, "%s type=%s\n", methodName.c_str(), (brookContext?"Brook":"Reference") );
(void) fflush( log );
BrookPlatform platform( 32, "cal", log );
Platform* platform;
if( brookContext ){
platform = new BrookPlatform( 32, "cal", log );
} else {
platform = new ReferencePlatform();
}
System* system = new System( numberOfAtoms, 0);
system->setAtomMass(0, 2.0);
system->setAtomMass(1, 2.0);
// double temperature, double frictionCoeff, double stepSize
LangevinIntegrator* integrator = new LangevinIntegrator(0, 0.1, 0.01);
*outIntegrator = integrator;
System system(numberOfAtoms, 0);
system.setAtomMass(0, 2.0);
system.setAtomMass(1, 2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 1, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1.5, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
system->addForce(forceField);
OpenMMContext* context = new OpenMMContext( *system, *integrator, *platform );
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
context->setPositions(positions);
return context;
}
void testLangevinSingleBond() {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "LangevinSingleBond";
static const int debug = 1;
FILE* log = stdout;
// ---------------------------------------------------------------------------------------
(void) fprintf( log, "%s\n", methodName.c_str() );
(void) fflush( log );
LangevinIntegrator* langevinIntegrator;
OpenMMContext* context = testLangevinSingleBondSetup( 1, &langevinIntegrator );
// This is simply a damped harmonic oscillator, so compare it to the analytical solution.
double freq = std::sqrt(1-0.05*0.05);
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Velocities);
double freq = std::sqrt(1-0.05*0.05);
int numberOfIterations = 1000;
for (int i = 0; i < numberOfIterations; ++i) {
State state = context->getState( State::Positions | State::Velocities );
double time = state.getTime();
double expectedDist = 1.5+0.5*std::exp(-0.05*time)*std::cos(freq*time);
// printf( "%s %d time=%.5e pos=[%.5e %.5e]\n", methodName.c_str(), i, time, -0.5*expectedDist, state.getPositions()[0] );
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedDist, 0, 0), state.getPositions()[0], 0.02);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedDist, 0, 0), state.getPositions()[1], 0.02);
double expectedSpeed = -0.5*std::exp(-0.05*time)*(0.05*std::cos(freq*time)+freq*std::sin(freq*time));
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedSpeed, 0, 0), state.getVelocities()[0], 0.02);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedSpeed, 0, 0), state.getVelocities()[1], 0.02);
integrator.step(1);
langevinIntegrator->step(1);
}
// Not set the friction to a tiny value and see if it conserves energy.
/*
integrator.setFriction(5e-5);
context.setPositions(positions);
State state = context.getState(State::Energy);
......@@ -1735,6 +2147,10 @@ void testLangevinSingleBond() {
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
*/
delete langevinIntegrator;
delete context;
}
void testLangevinTemperature() {
......@@ -1855,7 +2271,6 @@ int main( ){
// testBrookNonBonded( );
/*
testBrookStreams( 50 );
testBrookStreams( 63 );
testBrookStreams( 64 );
......@@ -1863,14 +2278,18 @@ int main( ){
testBrookBonds();
testBrookAngles();
testBrookRBTorsions();
testBrookPeriodicTorsions();
// testBrookRBTorsions();
testBrookCoulomb();
testBrookLJ();
testBrookExclusionsAnd14();
*/
testLangevinSingleBond();
testObcForce();
testObcSingleAtom();
testObcEConsistentForce();
//testForce(platform);
} catch( const exception& e ){
cout << "exception: " << e.what() << endl;
......
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