Skip to content

GitLab

Commit 6ccf6169 authored by Peter Eastman's avatar Peter Eastman
Browse files

Removed Brook platform

parent 82e0bd2f
Hide whitespace changes
Inline Side-by-side
Showing with 0 additions and 1958 deletions
platforms/brook/tests/TestBrookLangevinIntegrator.cpp deleted 100644 → 0
View file @ 82e0bd2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include <vector>
#include "../../../tests/AssertionUtilities.h"
#include "BrookPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/Context.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonBondedForce.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "../src/sfmt/SFMT.h"
#include "../../reference/src/SimTKUtilities/SimTKOpenMMRealType.h"
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
static Context* testLangevinSingleBondSetup( int brookContext, LangevinIntegrator** outIntegrator, FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "LangevinSingleBondSetup";
int PrintOn = 1;
int numberOfParticles = 2;
double mass = 2.0;
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
if( PrintOn ){
(void) fprintf( log, "%s type=%s\n", methodName.c_str(), (brookContext?"Brook":"Reference") );
(void) fflush( log );
}
Platform* platform;
if( brookContext ){
platform = new BrookPlatform( 32, "cal", log );
//platform = new BrookPlatform( 32, "cpu", log );
} else {
platform = new ReferencePlatform();
}
System* system = new System;
system->addParticle(mass);
system->addParticle(mass);
// double temperature, double frictionCoeff, double stepSize
LangevinIntegrator* integrator = new LangevinIntegrator(0, 0.1, 0.001);
integrator->setConstraintTolerance(1e-5);
*outIntegrator = integrator;
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system->addForce(forceField);
Context* context = new Context( *system, *integrator, *platform );
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
context->setPositions(positions);
return context;
}
void testLangevinSingleBond( FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "LangevinSingleBond";
int PrintOn = 1;
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
if( PrintOn ){
(void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
}
LangevinIntegrator* langevinIntegrator;
Context* context = testLangevinSingleBondSetup( 1, &langevinIntegrator, log );
// This is simply a damped harmonic oscillator, so compare it to the analytical solution.
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);
Vec3 pos1 = state.getPositions()[0];
Vec3 pos2 = state.getPositions()[1];
if( PrintOn > 1 ){
(void) fprintf( log, "%s %d time=%.5e expD=%.5e pos=[%.5f %.5f %.5f] [%.5f %.5f %.5f] ", methodName.c_str(), i, time, -0.5*expectedDist, pos1[0], pos1[1], pos1[2], pos2[0], pos2[1], pos2[2] );
(void) fflush( log );
}
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);
Vec3 vel1 = state.getVelocities()[0];
Vec3 vel2 = state.getVelocities()[1];
if( PrintOn > 1 ){
(void) fprintf( log, "expVel=%.5e vel=[%.5f %.5f %.5f] [%.5f %.5f %.5f]\n", -0.5*expectedSpeed, vel1[0], vel1[1], vel1[2], vel2[0], vel2[1], vel2[2] );
(void) fflush( stdout );
}
langevinIntegrator->step(1);
}
if( PrintOn ){
(void) fprintf( log, "%s 1 ok\n", methodName.c_str() ); fflush( log );
}
// Not set the friction to a tiny value and see if it conserves energy.
langevinIntegrator->setFriction(5e-5);
State state = context->getState(State::Energy);
double potentialEnergy = state.getPotentialEnergy();
double kineticEnergy = state.getKineticEnergy();
double initialEnergy = potentialEnergy + kineticEnergy;
if( PrintOn ){
(void) fprintf( log, "%s 2: initial energy: pot=%.5e ke=%.5e tot=%.5e\n", methodName.c_str(), potentialEnergy, kineticEnergy, initialEnergy );
(void) fflush( log );
}
for (int i = 0; i < 1000; ++i) {
state = context->getState(State::Energy);
double potentialEnergy = state.getPotentialEnergy();
double kineticEnergy = state.getKineticEnergy();
double energy = potentialEnergy + kineticEnergy;
if( PrintOn > 1 ){
(void) fprintf( log, "%s 2: energy: %d %.5e %.5e\n", methodName.c_str(), i, initialEnergy, energy, potentialEnergy, kineticEnergy );
(void) fflush( log );
}
ASSERT_EQUAL_TOL( initialEnergy, energy, 0.01);
langevinIntegrator->step(1);
}
if( PrintOn ){
(void) fprintf( log, "%s 2 ok\n", methodName.c_str() ); fflush( log );
}
//delete langevinIntegrator;
//delete context;
}
void testLangevinTemperature( FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "LangevinTemperature";
int PrintOn = 1;
const int numberOfParticles = 8;
double mass = 2.0;
const double temp = 100.0;
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
if( PrintOn ){
(void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
}
BrookPlatform platform( 32, "cal", log );
//ReferencePlatform platform;
System system;
LangevinIntegrator integrator(temp, 0.2, 0.002);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numberOfParticles; ++i){
system.addParticle(mass);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
CMMotionRemover* remover = new CMMotionRemover( 10 );
system.addForce(remover);
Context context(system, integrator, platform);
vector<Vec3> positions(numberOfParticles);
for (int i = 0; i < numberOfParticles; ++i){
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
}
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
int steps = 1000;
for( int i = 0; i < steps; ++i ){
State state = context.getState(State::Positions | State::Velocities | State::Energy);
//State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
if( PrintOn > 1 ){
(void) fprintf( log, "%s %d KE=%12.5e ttl=%12.5e\n",
methodName.c_str(), i, state.getKineticEnergy(), ke );
vector<Vec3> positions = state.getPositions();
vector<Vec3> velocities = state.getVelocities();
double com[3] = { 0.0, 0.0, 0.0 };
for( int ii = 0; ii < numberOfParticles; ii++ ){
com[0] += velocities[ii][0];
com[1] += velocities[ii][1];
com[2] += velocities[ii][2];
(void) fprintf( log, " %d q[%12.5e %12.5e %12.5e] v[%12.5e %12.5e %12.5e]\n", ii,
positions[ii][0], positions[ii][1], positions[ii][2],
velocities[ii][0], velocities[ii][1], velocities[ii][2] );
}
(void) fprintf( log, "VelCom[%12.5e %12.5e %12.5e]\n", com[0], com[1], com[2] );
(void)fflush( log );
}
integrator.step(1);
}
ke /= (double) steps;
double expected = 0.5*numberOfParticles*3.0*BOLTZ*temp;
double tol = 3*expected/std::sqrt(1000.0);
double diff = std::fabs( expected - ke );
if( PrintOn ){
(void) fprintf( log, "%s expected=%12.5e found=%12.5e diff=%12.5e tol=%12.5e\n", methodName.c_str(), expected, ke, diff, tol ); fflush( log );
}
ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
if( PrintOn ){
(void) fprintf( log, "%s ok\n", methodName.c_str(), expected, ke, diff, tol ); fflush( log );
}
/*
/tests/AssertionUtilities.h
#define ASSERT_EQUAL_TOL(expected, found, tol){
double _scale_ = std::fabs(expected) > 1.0 ? std::fabs(expected) : 1.0;
if (std::fabs((expected)-(found))/_scale_ > (t ol)) {std::stringstream details; details << "Expected "<<(expected)<<", found "<<(found); throwException(__FILE__, __LINE__, details.str());}};
ASSERT_EQUAL_TOL(expected, ke, tol );
*/
}
void testLangevinConstraints( FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "LangevinConstraints";
int PrintOn = 1;
double mass = 1.0;
// ---------------------------------------------------------------------------------------
if( PrintOn ){
(void) fprintf( log, "%s\n", methodName.c_str() );
(void) fflush( log );
}
BrookPlatform platform( 32, "cal", log );
const int numParticles = 8;
const int numConstraints = 4;
const double temp = 100.0;
// ReferencePlatform platform;
System system;
LangevinIntegrator integrator( temp, 2.0, 0.001 );
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(mass);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numConstraints; ++i){
system.addConstraint(2*i, 2*i+1, 1.0);
}
system.addForce(forceField);
CMMotionRemover* remover = new CMMotionRemover();
system.addForce(remover);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
init_gen_rand(0);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions);
for (int j = 0; j < numConstraints; ++j) {
int particle1, particle2;
double distance;
system.getConstraintParameters(j, particle1, particle2, distance);
Vec3 p1 = state.getPositions()[particle1];
Vec3 p2 = state.getPositions()[particle2];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
if( PrintOn > 1 ){
(void) fprintf( log, "%s %d %d dist=%12.5e %12.5e ok\n", methodName.c_str(), i, j, dist, fabs(dist-1.0) ); fflush( log );
}
ASSERT_EQUAL_TOL(1.0, dist, 2e-3);
}
integrator.step(1);
}
if( PrintOn ){
(void) fprintf( log, "%s ok\n", methodName.c_str() ); fflush( log );
}
}
int main( ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testBrookLangevinIntegrator";
FILE* log = stdout;
// ---------------------------------------------------------------------------------------
(void) fflush( stdout );
(void) fflush( stderr );
try {
// problem w/ testLangevinTemperature: T's not stable; however appears ok for bigger systems
// testLangevinTemperature( log );
testLangevinSingleBond( log );
testLangevinConstraints( log );
} catch( const exception& e ){
(void) fprintf( log, "Exception %s %.s\n", methodName.c_str(), e.what() ); (void) fflush( log );
return 1;
}
(void) fprintf( log, "\n%s done\n", methodName.c_str() ); (void) fflush( log );
return 0;
}
platforms/brook/tests/TestBrookNonbondedForce.cpp deleted 100644 → 0
View file @ 82e0bd2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
/**
* This tests the Brook harmonic angle bond force/energy
*/
#include <vector>
#include "../../../tests/AssertionUtilities.h"
#include "BrookPlatform.h"
#include "openmm/Context.h"
#include "openmm/NonbondedForce.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#define PI_M 3.141592653589
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testBrookCoulomb( FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "Coulomb";
int PrintOn = 0;
int numberOfParticles = 2;
double mass = 2.0;
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
if( PrintOn ){
(void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
}
BrookPlatform platform( 32, "cal", log );
System system;
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator( 0, 0.1, 0.01 );
// int index, double charge, double radius, double depth
NonbondedForce* forceField = new NonbondedForce();
forceField->addParticle(0.5, 1, 0);
forceField->addParticle(-1.5, 1, 0);
system.addForce(forceField);
//(void) fprintf( log, "%s: Calling context\n", methodName.c_str() );
//(void) fflush( log );
Context context(system, integrator, platform);
vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
context.setPositions(positions);
//(void) fprintf( log, "%s :Calling getState\n", methodName.c_str() );
//(void) fflush( log );
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
if( PrintOn ){
(void) fprintf( log, "\nCoulomb forces\n");
for( int ii = 0; ii < numberOfParticles; ii++ ){
(void) fprintf( log, "%d [%.5e %.5e %.5e]\n", ii, forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
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);
if( PrintOn ){
(void) fprintf( log, "Coulomb forces ok\n"); fflush( log );
}
// delete forceField;
}
void testBrookLJ( FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "LJ";
int PrintOn = 0;
int numberOfParticles = 2;
double mass = 2.0;
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
if( PrintOn ){
(void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
}
BrookPlatform platform( 32, "cal", log );
// ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator( 0, 0.1, 0.01 );
// int index, double charge, double radius, double depth
NonbondedForce* forceField = new NonbondedForce();
forceField->addParticle(0, 1.2, 1);
forceField->addParticle(0, 1.4, 2);
system.addForce(forceField);
//(void) fprintf( log, "%s: Calling context\n", methodName.c_str() );
//(void) fflush( log );
Context context(system, integrator, platform);
vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
context.setPositions(positions);
//(void) fprintf( log, "%s :Calling getState\n", methodName.c_str() );
//(void) fflush( log );
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
if( PrintOn ){
(void) fprintf( log, "LJ forces\n");
for( int ii = 0; ii < numberOfParticles; ii++ ){
(void) fprintf( log, "%d [%.5e %.5e %.5e]\n", ii, forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
double x = 1.3/2.0;
double eps = sqrt( 2.0 );
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);
if( PrintOn ){
(void) fprintf( log, "LJ forces ok\n"); fflush( log );
}
}
void testBrookExclusionsAnd14( FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "ExclusionsAnd14";
int numberOfParticles = 5;
int PrintOn = 0;
double mass = 2.0;
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
if( PrintOn ){
(void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
}
BrookPlatform platform( 32, "cpu", log );
//ReferencePlatform platform;
System system;
LangevinIntegrator integrator( 0, 0.1, 0.01 );
// int index, double charge, double radius, double depth
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < numberOfParticles; i++) {
system.addParticle(1.0);
nonbonded->addParticle(0, 1.5, 0);
}
vector<pair<int, int> > bonds;
bonds.push_back(pair<int, int>(0, 1));
bonds.push_back(pair<int, int>(1, 2));
bonds.push_back(pair<int, int>(2, 3));
bonds.push_back(pair<int, int>(3, 4));
nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
int first14, second14;
for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon;
nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
first14 = i;
if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
second14 = i;
}
system.addForce(nonbonded);
//(void) fprintf( log, "%s: Calling context\n", methodName.c_str() );
//(void) fflush( log );
Context context(system, integrator, platform);
vector<Vec3> positions(numberOfParticles);
const double r = 1.0;
positions[0] = Vec3(0, 0, 0);
for( int ii = 1; ii < numberOfParticles; ii++ ){
positions[ii] = Vec3(r, 0, 0);
}
for( int ii = 1; ii < numberOfParticles; ii++ ){
// Test LJ forces
vector<Vec3> positions(5);
const double r = 1.0;
for (int j = 0; j < 5; ++j) {
nonbonded->setParticleParameters(j, 0, 1.5, 0);
positions[j] = Vec3(0, j, 0);
}
nonbonded->setParticleParameters(0, 0, 1.5, 1);
nonbonded->setParticleParameters(ii, 0, 1.5, 1);
nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, ii == 3 ? 0.5 : 0.0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
positions[ii] = Vec3(r, 0, 0);
context.reinitialize();
context.setPositions(positions);
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
double x = 1.5/r;
double eps = 1.0;
double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
double energy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
if( ii == 3 ){
force *= 0.5;
energy *= 0.5;
}
if( ii < 3 ){
force = 0;
energy = 0;
}
if( PrintOn ){
(void) fprintf( log, "14 LJ forces ii=%d F=%.6e\n", ii, force );
for( int jj = 0; jj < numberOfParticles; jj++ ){
(void) fprintf( log, "%d [%.5e %.5e %.5e]\n", jj, forces[jj][0], forces[jj][1], forces[jj][2] );
}
(void) fflush( log );
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[ii], TOL);
if( PrintOn ){
(void) fprintf( log, "14 LJ forces ok for index=%d\n\n", ii );
(void) fflush( log );
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
// Test Coulomb forces
nonbonded->setParticleParameters( 0, 2, 1.5, 0 );
nonbonded->setParticleParameters( ii, 2, 1.5, 0 );
nonbonded->setExceptionParameters( first14, 0, 3, ii == 3 ? 4/1.2 : 0, 1.5, 0 );
nonbonded->setExceptionParameters( second14, 1, 4, 0, 1.5, 0 );
context.reinitialize();
context.setPositions(positions);
state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces2 = state.getForces();
force = 138.935485*4/(r*r);
energy = 138.935485*4/r;
if( ii == 3 ){
force /= 1.2;
energy /= 1.2;
}
if( ii < 3 ){
force = 0;
energy = 0;
}
if( PrintOn ){
(void) fprintf( log, "14 Coulomb forces ii=%d F=%.6e\n", ii, force );
for( int jj = 0; jj < numberOfParticles; jj++ ){
(void) fprintf( log, "%d [%.5e %.5e %.5e]\n", jj, forces2[jj][0], forces2[jj][1], forces2[jj][2] );
}
(void) fflush( log );
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[ii], TOL);
if( PrintOn ){
(void) fprintf( log, "14 Coulomb forces ok for index=%d\n\n", ii );
(void) fflush( log );
}
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
if( PrintOn ){
(void) fprintf( log, "ExclusionsAnd14 ok\n"); fflush( log );
}
}
int main( ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testBrookNonBonded";
FILE* log = stdout;
// ---------------------------------------------------------------------------------------
(void) fflush( stdout );
(void) fflush( stderr );
try {
testBrookCoulomb( log );
testBrookLJ( log );
testBrookExclusionsAnd14( log );
} catch( const exception& e ){
(void) fprintf( log, "Exception %s %.s\n", methodName.c_str(), e.what() ); (void) fflush( log );
return 1;
}
(void) fprintf( log, "\n%s done\n", methodName.c_str() ); (void) fflush( log );
return 0;
}
platforms/brook/tests/TestBrookPeriodicTorsionForce.cpp deleted 100644 → 0
View file @ 82e0bd2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
/**
* This tests the Brook periodic torsion bond force/energy
*/
#include "../../../tests/AssertionUtilities.h"
#include "BrookPlatform.h"
#include "openmm/Context.h"
#include "openmm/PeriodicTorsionForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include <vector>
#define PI_M 3.141592653589
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testBrookPeriodicTorsions( FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "PeriodicTorsions";
int PrintOn = 0;
int numberOfParticles = 4;
double mass = 2.0;
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
if( PrintOn ){
(void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
}
BrookPlatform platform( 32, "cal", log );
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator( 0, 0.1, 0.01 );
// int numParticles, int numBonds, int numAngles, int numPeriodicTorsions, int numRBTorsions
PeriodicTorsionForce* forceField = new PeriodicTorsionForce();
// int index, int atom1, int atom2, int atom3, double angle, double k
forceField->addTorsion(0, 1, 2, 3, 2, PI_M/3, 1.1);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
positions[3] = Vec3(1, 0, 2);
context.setPositions(positions);
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
if( PrintOn ){
(void) fprintf( log, "Periodic torsion bond forces\n");
for( int ii = 0; ii < numberOfParticles; ii++ ){
(void) fprintf( log, "%d [%.5e %.5e %.5e]\n", ii, forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
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), tolerance );
ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), tolerance );
if( PrintOn ){
(void) fprintf( log, "Periodic torsion bond forces ok -- tolerance=%.2e\n", tolerance); (void) fflush( log );
}
}
int main( ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testBrookPeriodicTorsions";
FILE* log = stdout;
// ---------------------------------------------------------------------------------------
(void) fflush( stdout );
(void) fflush( stderr );
try {
testBrookPeriodicTorsions( log );
} catch( const exception& e ){
(void) fprintf( log, "Exception %s %.s\n", methodName.c_str(), e.what() ); (void) fflush( log );
return 1;
}
(void) fprintf( log, "\n%s done\n", methodName.c_str() ); (void) fflush( log );
return 0;
}
platforms/brook/tests/TestBrookRBTorsionForce.cpp deleted 100644 → 0
View file @ 82e0bd2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
/**
* This tests the Brook RB torsion bond force/energy
*/
#include "../../../tests/AssertionUtilities.h"
#include "BrookPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/Context.h"
#include "openmm/RBTorsionForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include <vector>
#define PI_M 3.141592653589
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testBrookRBTorsions( FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "RBTorsions";
int PrintOn = 0;
int numberOfParticles = 4;
double mass = 2.0;
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
if( PrintOn ){
(void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
}
BrookPlatform platform( 32, "cal", log );
//ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator( 0, 0.1, 0.01 );
RBTorsionForce* forceField = new RBTorsionForce();
forceField->addTorsion(0, 1, 2, 3, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(1, 0, 0);
positions[3] = Vec3(1, 1, 1);
context.setPositions(positions);
State state = context.getState( State::Forces | State::Energy );
const vector<Vec3>& forces = state.getForces();
if( PrintOn ){
(void) fprintf( log, "RB torsion bond forces\n");
for( int ii = 0; ii < numberOfParticles; ii++ ){
(void) fprintf( log, "%d [%.5e %.5e %.5e]\n", ii, forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
double psi = 0.25*PI_M - PI_M;
double torque = 0.0;
for (int i = 1; i < 6; ++i) {
double c = 0.1*(i+1);
torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
}
if( PrintOn ){
(void) fprintf( log, "RB torsion bond expected forces\n");
(void) fprintf( log, "0 [0.0 0.0 %.5e]\n", torque );
(void) fprintf( log, "3 [0.0 %.5e %.5e]\n", 0.5*torque, -0.5*torque );
(void) fflush( log );
}
double tolerance = 0.001;
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), 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 );
if( PrintOn ){
(void) fprintf( log, "RB torsion bond forces ok tolerance=%.2e\n", tolerance); fflush( log );
}
}
int main( ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testBrookRbTorsion";
FILE* log = stdout;
// ---------------------------------------------------------------------------------------
(void) fflush( stdout );
(void) fflush( stderr );
try {
testBrookRBTorsions( log );
} catch( const exception& e ){
(void) fprintf( log, "Exception %s %.s\n", methodName.c_str(), e.what() ); (void) fflush( log );
return 1;
}
(void) fprintf( log, "\n%s done\n", methodName.c_str() ); (void) fflush( log );
return 0;
}
platforms/brook/tests/TestBrookStream.cpp deleted 100644 → 0
View file @ 82e0bd2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
/**
* This tests the Brook stream.
*/
#include "../../../tests/AssertionUtilities.h"
#include "BrookPlatform.h"
#include "BrookStreamFactory.h"
#include <iostream>
#include <vector>
#include <fstream>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testWriteRead( void ){
static const int ArraySz = 3000;
BrookPlatform platform;
// create and initialize arrays
float* array = new float[ArraySz];
float* saveArray = new float[ArraySz];
for( int ii = 0; ii < ArraySz; ii++ ){
array[ii] = (float) ii;
}
memset( saveArray, 0, sizeof( float )*ArraySz );
// get factory & create stream
const BrookStreamFactory& brookStreamFactory = dynamic_cast<const BrookStreamFactory&> (platform.getDefaultStreamFactory());
StreamImpl* testStream = brookStreamFactory.createStreamImpl( "TestStream", ArraySz, Stream::Float, platform );
// load & retreive data
testStream->loadFromArray( array );
testStream->saveToArray( saveArray );
// test for equality
for( int ii = 0; ii < ArraySz; ii++ ){
ASSERT_EQUAL( array[ii], saveArray[ii] );
}
delete[] saveArray;
delete[] array;
}
int main( ){
(void) fflush( stdout );
(void) fflush( stderr );
try {
testWriteRead();
} catch( const exception& e ){
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/brook/tests/TestBrookVerletIntegrator.cpp deleted 100644 → 0
View file @ 82e0bd2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
/**
* This tests the Brook harmonic angle bond force/energy
*/
#include <vector>
#include "../../../tests/AssertionUtilities.h"
#include "BrookPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/Context.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "../src/sfmt/SFMT.h"
#define PI_M 3.141592653589
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testVerletSingleBond( FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testVerletSingleBond";
int PrintOn = 1;
int numberOfParticles = 2;
double mass = 2.0;
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
if( PrintOn ){
(void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
}
BrookPlatform platform( 32, "cal", log );
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VerletIntegrator integrator(0.001);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
// CMMotionRemover* remover = new CMMotionRemover();
// system.addForce(remover);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
// This is simply a harmonic oscillator, so compare it to the analytical solution.
const double freq = 1.0;;
State state = context.getState(State::Energy);
const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
if( PrintOn ){
(void) fprintf( log, "%s Energy initialEnergy=%12.5e KE=%12.5e PE=%12.5e\n",
methodName.c_str(), initialEnergy,
state.getKineticEnergy(), state.getPotentialEnergy() );
(void) fflush( log );
}
for (int i = 0; i < 1000; ++i) {
state = context.getState(State::Positions | State::Velocities | State::Energy);
double time = state.getTime();
double expectedDist = 1.5+0.5*std::cos(freq*time);
Vec3 position0 = state.getPositions()[0];
Vec3 position1 = state.getPositions()[1];
if( PrintOn > 1 ){
(void) fprintf( log, "%s %d Pos expected=[%12.5e 0 0] actual=[%12.5e %12.5e %12.5e] [%12.5e %12.5e %12.5e]\n",
methodName.c_str(), i, -0.5*expectedDist,
position0[0], position0[1], position0[2],
position1[0], position1[1], position1[2] );
(void) fflush( log );
}
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*freq*std::sin(freq*time);
Vec3 velocity0 = state.getVelocities()[0];
Vec3 velocity1 = state.getVelocities()[1];
if( PrintOn > 1 ){
(void) fprintf( log, "%s %d Vel expected=[%12.5e 0 0] actual=[%12.5e %12.5e %12.5e] [%12.5e %12.5e %12.5e]\n",
methodName.c_str(), i, -0.5*expectedSpeed,
velocity0[0], velocity0[1], velocity0[2],
velocity1[0], velocity1[1], velocity1[2] );
(void) fflush( log );
}
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);
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if( PrintOn > 1 ){
(void) fprintf( log, "%s %d Energy initialEnergy=%12.5e actual=%12.5e KE=%12.5e PE=%12.5e\n",
methodName.c_str(), i, initialEnergy, energy,
state.getKineticEnergy(), state.getPotentialEnergy() );
(void) fflush( log );
}
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
if( PrintOn ){
(void) fprintf( log, "%s ok\n", methodName.c_str() );
(void) fflush( log );
}
}
void testVerletConstraints( FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testVerletConstraints";
int PrintOn = 1;
const int numParticles = 8;
const int numConstraints = numParticles/2;
double mass = 10.0;
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
if( PrintOn ){
(void) fprintf( log, "%s\n", methodName.c_str() ); (void) fflush( log );
}
//ReferencePlatform platform;
BrookPlatform platform( 32, "cal", log );
System system;
VerletIntegrator integrator(0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(mass);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numConstraints; ++i){
system.addConstraint(2*i, 2*i+1, 1.0);
}
system.addForce(forceField);
//CMMotionRemover* remover = new CMMotionRemover();
//system.addForce(remover);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
init_gen_rand(0);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
double tolerance = 0.002;
double maxDiff = -1.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy);
for (int j = 0; j < numConstraints; ++j) {
int particle1, particle2;
double distance;
system.getConstraintParameters(j, particle1, particle2, distance);
Vec3 p1 = state.getPositions()[particle1];
Vec3 p2 = state.getPositions()[particle2];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
double diff = fabs( distance - dist );
if( diff > maxDiff ){
maxDiff = diff;
}
if( PrintOn > 1 || diff > tolerance ){
(void) fprintf( log, "%s step=%d cnstrnt=%d p[%d %d] d=%.5e exptd=%.5e dif=%.5e [%.5e %.5e %.5e] [%.5e %.5e %.5e] mxDff=%.5e\n",
methodName.c_str(), i, j, particle1, particle2, dist, distance, diff,
p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], maxDiff); (void) fflush( log );
}
ASSERT_EQUAL_TOL(distance, dist, tolerance );
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if( PrintOn > 1 ){
(void) fprintf( log, "%s %d e[%.5e %.5e] ke=%.5e pe=%.5e\n",
methodName.c_str(), i, initialEnergy, energy, state.getKineticEnergy(), state.getPotentialEnergy() ); (void) fflush( log );
}
if( i == 1 ){
initialEnergy = energy;
} else if( i > 1 ){
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.5);
}
integrator.step(1);
}
if( PrintOn ){
(void) fprintf( log, "%s ok maxShakeDiff=%.5e tolerance=%.5e\n", methodName.c_str(), maxDiff, tolerance );
(void) fflush( log );
}
}
int main( ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testBrookVerletIntegrator";
FILE* log = stdout;
// ---------------------------------------------------------------------------------------
(void) fflush( stdout );
(void) fflush( stderr );
try {
testVerletSingleBond( log );
testVerletConstraints( log );
} catch( const exception& e ){
(void) fprintf( log, "Exception %s %.s\n", methodName.c_str(), e.what() ); (void) fflush( log );
return 1;
}
(void) fprintf( log, "\n%s done\n", methodName.c_str() ); (void) fflush( log );
return 0;
}
platforms/brook/tests/TstBrookGBSAOBCForce.cpp deleted 100644 → 0
View file @ 82e0bd2f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2009 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
/**
* This tests the Brook OBC force/energy
*/
#include <vector>
#include <fstream>
#include <iostream>
#include "../../../tests/AssertionUtilities.h"
#include "BrookPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/Context.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "../src/sfmt/SFMT.h"
#include "../../reference/src/SimTKUtilities/SimTKOpenMMRealType.h"
typedef std::vector<std::string> StringVector;
typedef StringVector::iterator StringVectorI;
typedef StringVector::const_iterator StringVectorCI;
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
/**---------------------------------------------------------------------------------------
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
--------------------------------------------------------------------------------------- */
char* strsepLocal( 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 );
}
}
/**---------------------------------------------------------------------------------------
Tokenize a string (static method) (Simbios)
@param lineBuffer string to tokenize
@param tokenArray upon return vector of tokens
@param delimiter token delimter
@return number of args
--------------------------------------------------------------------------------------- */
int tokenizeString( char* lineBuffer, StringVector& tokenArray, const std::string delimiter ){
// ---------------------------------------------------------------------------------------
// static const std::string methodName = "\nSimTKOpenMMUtilities::tokenizeString";
// ---------------------------------------------------------------------------------------
// (void) fprintf( stdout, "\nIn SimTKOpenMMUtilities::tokenizeString <%s>", lineBuffer );
// (void) fflush( stdout );
char *ptr_c = NULL;
for( ; (ptr_c = strsepLocal( &lineBuffer, delimiter.c_str() )) != NULL; ){
if( *ptr_c ){
tokenArray.push_back( std::string( ptr_c ) );
}
}
return (int) tokenArray.size();
}
static Context* testObcForceSetup( int numParticles, int brookContext, FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testObcForceSetup";
// ---------------------------------------------------------------------------------------
Platform* platform;
if( brookContext ){
platform = new BrookPlatform( 32, "cal", log );
} else {
platform = new ReferencePlatform();
}
System* system = new System();
for (int i = 0; i < numParticles; i++)
system->addParticle(1.0);
LangevinIntegrator* integrator = new LangevinIntegrator(0, 0.1, 0.01);
GBSAOBCForce* forceField = new GBSAOBCForce();
for (int i = 0; i < numParticles; ++i){
// charge radius scalingFactor
forceField->addParticle(i%2 == 0 ? -1 : 1, 0.15, 1);
//forceField->setParticleParameters(i, i%2 == 0 ? -1 : 1, 1.5, 1);
}
system->addForce(forceField);
Context* context = new Context( *system, *integrator, *platform );
// Set random positions for all the atoms.
vector<Vec3> positions(numParticles);
init_gen_rand(0);
for (int i = 0; i < numParticles; ++i){
positions[i] = Vec3(1.0*genrand_real2(), 1.0*genrand_real2(), 1.0*genrand_real2());
}
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 Context* testObcForceFileSetup( std::string fileName, int brookContext, int* numParticles, FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testObcForceFileSetup";
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( log ){
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 && log ){
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 numberOfParticles = atoi( (*ii).c_str() );
if( log ){
(void) fprintf( log, "\n%s %d\n", methodName.c_str(), numberOfParticles );
}
*numParticles = numberOfParticles;
lineCount++;
System* system = new System();
for (int i = 0; i < numberOfParticles; i++)
system->addParticle(1.0);
GBSAOBCForce* forceField = new GBSAOBCForce();
vector<Vec3> positions(numberOfParticles);
int index = 0;
for (int i = 0; i < numberOfParticles; ++i){
(void) fgets( buffer, bufferSize, readFile );
StringVector tokens;
tokenizeString( buffer, tokens, delimiter );
if( tokens.size() < 8 && log ){
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 = (float) atof( (*ii).c_str() ); ii++;
float coordY = (float) atof( (*ii).c_str() ); ii++;
float coordZ = (float) atof( (*ii).c_str() ); ii++;
float radius = (float) atof( (*ii).c_str() ); ii++;
float scalingFactor = (float) atof( (*ii).c_str() ); ii++;
float charge = (float) atof( (*ii).c_str() ); ii++;
float bornRadi = (float) atof( (*ii).c_str() ); ii++;
positions[index++] = Vec3( coordX, coordY, coordZ );
// charge radius scalingFactor
forceField->addParticle( charge, radius, scalingFactor );
if( log ){
(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++;
}
if( log ){
(void) fflush( log );
}
system->addForce(forceField);
Context* context = new Context( *system, *integrator, *platform );
// Set positions for all the atoms.
context->setPositions(positions);
return context;
}
void testObcForce( FILE* log, char* testInputFileName ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testObcForce";
static const int maxPrint = 10;
static int PrintOn = 0;
float epsilon = 1.0e-04f;
int totalErrors = 0;
int numParticles = 10;
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
log = PrintOn ? log : NULL;
//Context* context = testObcForceSetup( numParticles, 0, log );
//Context* brookContext = testObcForceSetup( numParticles, 1, log );
Context* context = testObcForceFileSetup( std::string( testInputFileName ), 0, &numParticles, log );
//Context* context = NULL;
Context* brookContext = testObcForceFileSetup( std::string( testInputFileName ), 1, &numParticles, log );
//Context* brookContext = NULL;
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();
}
if( PrintOn > 1 ){
(void) fprintf( log, "%s OBC forces [%s %s]\n", methodName.c_str(), (context?"Ref":""), (brookContext?"Brook":"") );
for( int ii = 0; ii < numParticles; 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 < numParticles; ++i) {
Vec3 f = forces[i];
Vec3 fBrook = brookForces[i];
ASSERT_EQUAL_VEC( f, fBrook, tolerance );
}
}
if( PrintOn ){
(void) fprintf( log, "testObcForce ok w/ tolerance=%.3e\n", tolerance );
(void) fflush( log );
}
}
void testObcSingleParticle( FILE* log ){
// ---------------------------------------------------------------------------------------
int PrintOn = 0;
static const std::string methodName = "testObcSingleParticle";
// ---------------------------------------------------------------------------------------
PrintOn = log ? PrintOn : 0;
BrookPlatform platform( 32, "cal", log );
//ReferencePlatform platform;
System system;
system.addParticle(2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* forceField = new GBSAOBCForce();
forceField->addParticle(0.5, 0.15, 1);
system.addForce(forceField);
Context 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);
if( PrintOn ){
(void) fprintf( log, "%s ok\n", methodName.c_str() );
(void) fflush( log );
}
}
void testObcEConsistentForce( FILE* log ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testObcEConsistentForce";
static int PrintOn = 0;
// ---------------------------------------------------------------------------------------
BrookPlatform platform( 32, "cal", log );
//ReferencePlatform platform;
const int numParticles = 10;
System system;
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* forceField = new GBSAOBCForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
forceField->addParticle(i%2 == 0 ? -1 : 1, 0.15, 1);
}
system.addForce(forceField);
Context context(system, integrator, platform);
// Set random positions for all the atoms.
vector<Vec3> positions(numParticles);
init_gen_rand(0);
for (int i = 0; i < numParticles; ++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 < numParticles; ++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 < numParticles; ++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)
if( PrintOn ){
(void) fprintf( log, "%s ok %.8e %.8e %.8e %.8e\n", methodName.c_str(),
state2.getPotentialEnergy(), state.getPotentialEnergy(), norm,
(state2.getPotentialEnergy()-state.getPotentialEnergy())/delta ); (void) fflush( log );
}
}
int main( int argc, char* argv[] ){
// ---------------------------------------------------------------------------------------
static const std::string methodName = "testBrookObcGbsa";
FILE* log = stdout;
// ---------------------------------------------------------------------------------------
(void) fflush( stdout ); (void) fflush( stderr );
try {
// skip test if input file name is not provided
// but let user know
if( argc > 1 ){
testObcForce( log, argv[1] );
} else {
(void) fprintf( log, "%s missing input file (=ObcInfo.txt?): skipping testObcForce() argc=%d\n", methodName.c_str(), argc );
for( int ii = 0; ii < argc; ii++ ){
(void) fprintf( log, "%d <%s>\n", ii, argv[ii] );
}
(void) fflush( log );
}
testObcSingleParticle( log );
testObcEConsistentForce( log );
} catch( const exception& e ){
(void) fprintf( log, "Exception %s %.s\n", methodName.c_str(), e.what() ); (void) fflush( log );
return 1;
}
(void) fprintf( log, "\n%s done\n", methodName.c_str() ); (void) fflush( log );
return 0;
}
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