Skip to content

GitLab

"plugins/vscode:/vscode.git/clone" did not exist on "21ab511a59dc28926ebc1c335ad134b79f6a0a7d"
Commit 8b5d3dc0 authored by Peter Eastman's avatar Peter Eastman
Browse files

Finished refactoring test cases

parent 2f553a66
Hide whitespace changes
Inline Side-by-side
Showing with 1418 additions and 3727 deletions
platforms/cuda/tests/TestCudaVariableLangevinIntegrator.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,311 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of VariableLangevinIntegrator.
*/
#include "CudaTests.h"
#include "TestVariableLangevinIntegrator.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VariableLangevinIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CudaPlatform platform;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VariableLangevinIntegrator integrator(0, 0.1, 1e-6);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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 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 time = state.getTime();
double expectedDist = 1.5+0.5*std::exp(-0.05*time)*std::cos(freq*time);
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);
}
// Now 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);
double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
for (int i = 0; i < 1000; ++i) {
state = context.getState(State::Energy);
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
integrator.step(1);
}
}
void testTemperature() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 5.0, 5e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; ++i)
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Let it equilibrate.
integrator.step(5000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 5000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(5);
}
ke /= 5000;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.1);
}
void testConstraints() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 2.0, 1e-5);
integrator.setConstraintTolerance(1e-5);
integrator.setRandomNumberSeed(0);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numParticles-1; ++i)
system.addConstraint(i, i+1, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-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 < numParticles-1; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
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]));
ASSERT_EQUAL_TOL(1.0, dist, 2e-5);
}
integrator.step(1);
}
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VariableLangevinIntegrator integrator(300.0, 2.0, 0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
void testRandomSeed() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 2.0, 1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(0, 0, 0);
}
// Try twice with the same random seed.
integrator.setRandomNumberSeed(5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state1 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state2 = context.getState(State::Positions);
// Try twice with a different random seed.
integrator.setRandomNumberSeed(10);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state3 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state4 = context.getState(State::Positions);
// Compare the results.
for (int i = 0; i < numParticles; i++) {
for (int j = 0; j < 3; j++) {
ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
}
}
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
int numParticles = 0;
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
++numParticles;
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableLangevinIntegrator integrator(temp, 6.0, 1e-4);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(2.0);
// Make sure the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 1000; ++i) {
double t = 2.0 + 0.02 * (i + 1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
}
ke /= 1000;
double expected = 1.5 * numParticles * BOLTZ * temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.01);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testSingleBond();
testTemperature();
testConstraints();
testConstrainedMasslessParticles();
testRandomSeed();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cuda/tests/TestCudaVariableVerletIntegrator.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,289 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of VariableVerletIntegrator.
*/
#include "CudaTests.h"
#include "TestVariableVerletIntegrator.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VariableVerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CudaPlatform platform;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VariableVerletIntegrator integrator(1e-6);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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();
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);
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);
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();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
integrator.step(1);
}
}
void testConstraints() {
const int numParticles = 8;
const int numConstraints = 5;
const double temp = 100.0;
System system;
VariableVerletIntegrator integrator(1e-5);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(1, 2, 1.0);
system.addConstraint(2, 3, 1.0);
system.addConstraint(4, 5, 1.0);
system.addConstraint(6, 7, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
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]));
ASSERT_EQUAL_TOL(distance, dist, 1e-4);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
double finalTime = context.getState(State::Positions).getTime();
ASSERT(finalTime > 0.1);
// Now try the stepTo() method.
finalTime += 0.5;
integrator.stepTo(finalTime);
ASSERT_EQUAL(finalTime, context.getState(State::Positions).getTime());
}
void testConstrainedClusters() {
const int numParticles = 7;
System system;
VariableVerletIntegrator integrator(1e-5);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i > 1 ? 1.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(0, 2, 1.0);
system.addConstraint(0, 3, 1.0);
system.addConstraint(0, 4, 1.0);
system.addConstraint(1, 5, 1.0);
system.addConstraint(1, 6, 1.0);
system.addConstraint(2, 3, sqrt(2.0));
system.addConstraint(2, 4, sqrt(2.0));
system.addConstraint(3, 4, sqrt(2.0));
system.addConstraint(5, 6, sqrt(2.0));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(-1, 0, 0);
positions[3] = Vec3(0, 1, 0);
positions[4] = Vec3(0, 0, 1);
positions[5] = Vec3(2, 0, 0);
positions[6] = Vec3(1, 1, 0);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i)
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
for (int j = 0; j < system.getNumConstraints(); ++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]));
ASSERT_EQUAL_TOL(distance, dist, 2e-5);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
ASSERT(context.getState(State::Positions).getTime() > 0.1);
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VariableVerletIntegrator integrator(0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableVerletIntegrator integrator(1e-5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(1.0);
// Simulate it and see whether energy remains constant.
State state0 = context.getState(State::Energy);
double initialEnergy = state0.getKineticEnergy() + state0.getPotentialEnergy();
for (int i = 0; i < 20; i++) {
double t = 1.0 + 0.05*(i+1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
double energy = state.getKineticEnergy() + state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
}
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testSingleBond();
testConstraints();
testConstrainedClusters();
testConstrainedMasslessParticles();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cuda/tests/TestCudaVerletIntegrator.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,225 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of VerletIntegrator.
*/
#include "CudaTests.h"
#include "TestVerletIntegrator.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CudaPlatform platform;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VerletIntegrator integrator(0.01);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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();
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);
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);
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();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
ASSERT_EQUAL_TOL(10.0, context.getState(0).getTime(), 1e-5);
}
void testConstraints() {
const int numParticles = 8;
const int numConstraints = 5;
const double temp = 100.0;
System system;
VerletIntegrator integrator(0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(1, 2, 1.0);
system.addConstraint(2, 3, 1.0);
system.addConstraint(4, 5, 1.0);
system.addConstraint(6, 7, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
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]));
ASSERT_EQUAL_TOL(distance, dist, 1e-4);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
void runPlatformTests() {
}
void testConstrainedClusters() {
const int numParticles = 7;
System system;
VerletIntegrator integrator(0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i > 1 ? 1.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(0, 2, 1.0);
system.addConstraint(0, 3, 1.0);
system.addConstraint(0, 4, 1.0);
system.addConstraint(1, 5, 1.0);
system.addConstraint(1, 6, 1.0);
system.addConstraint(2, 3, sqrt(2.0));
system.addConstraint(2, 4, sqrt(2.0));
system.addConstraint(3, 4, sqrt(2.0));
system.addConstraint(5, 6, sqrt(2.0));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(-1, 0, 0);
positions[3] = Vec3(0, 1, 0);
positions[4] = Vec3(0, 0, 1);
positions[5] = Vec3(2, 0, 0);
positions[6] = Vec3(1, 1, 0);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i)
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
for (int j = 0; j < system.getNumConstraints(); ++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]));
ASSERT_EQUAL_TOL(distance, dist, 2e-5);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VerletIntegrator integrator(0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testSingleBond();
testConstraints();
testConstrainedClusters();
testConstrainedMasslessParticles();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cuda/tests/TestCudaVirtualSites.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2012-2014 Stanford University and the Authors. *
* Portions copyright (c) 2012-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,394 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of virtual sites.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/CustomBondForce.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/VirtualSite.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CudaPlatform platform;
/**
* Check that massless particles are handled correctly.
*/
void testMasslessParticle() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
CustomBondForce* bonds = new CustomBondForce("-1/r");
system.addForce(bonds);
vector<double> params;
bonds->addBond(0, 1, params);
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
vector<Vec3> velocities(2);
velocities[0] = Vec3(0, 0, 0);
velocities[1] = Vec3(0, 1, 0);
context.setVelocities(velocities);
// The second particle should move in a circular orbit around the first one.
// Compare it to the analytical solution.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Velocities | State::Forces);
double time = state.getTime();
ASSERT_EQUAL_VEC(Vec3(0,0,0), state.getPositions()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(0,0,0), state.getVelocities()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(cos(time), sin(time), 0), state.getPositions()[1], 0.01);
ASSERT_EQUAL_VEC(Vec3(-sin(time), cos(time), 0), state.getVelocities()[1], 0.01);
integrator.step(1);
}
}
/**
* Test a TwoParticleAverageSite virtual site.
*/
void testTwoParticleAverage() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.8, 0.2));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.8+pos[1]*0.2, pos[2], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.1+0.3*0.8, 0, 0), state.getForces()[0], 1e-4);
ASSERT_EQUAL_VEC(Vec3(0.2+0.3*0.2, 0, 0), state.getForces()[1], 1e-4);
integrator.step(1);
}
}
/**
* Test a ThreeParticleAverageSite virtual site.
*/
void testThreeParticleAverage() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new ThreeParticleAverageSite(0, 1, 2, 0.2, 0.3, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.2+pos[1]*0.3+pos[2]*0.5, pos[3], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4*0.2, 0, 0), state.getForces()[0], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.2+0.4*0.3, 0, 0), state.getForces()[1], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.3+0.4*0.5, 0, 0), state.getForces()[2], 1e-5);
integrator.step(1);
}
}
/**
* Test an OutOfPlaneSite virtual site.
*/
void testOutOfPlane() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new OutOfPlaneSite(0, 1, 2, 0.3, 0.4, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 v12 = pos[1]-pos[0];
Vec3 v13 = pos[2]-pos[0];
Vec3 cross = v12.cross(v13);
ASSERT_EQUAL_VEC(pos[0]+v12*0.3+v13*0.4+cross*0.5, pos[3], 1e-5);
const vector<Vec3>& f = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0.1+0.2+0.3+0.4, 0, 0), f[0]+f[1]+f[2], 1e-5);
Vec3 f2(0.4*0.3, 0.4*0.5*v13[2], -0.4*0.5*v13[1]);
Vec3 f3(0.4*0.4, -0.4*0.5*v12[2], 0.4*0.5*v12[1]);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4, 0, 0)-f2-f3, f[0], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.2, 0, 0)+f2, f[1], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.3, 0, 0)+f3, f[2], 1e-5);
integrator.step(1);
}
}
/**
* Test a LocalCoordinatesSite virtual site.
*/
void testLocalCoordinates() {
const Vec3 originWeights(0.2, 0.3, 0.5);
const Vec3 xWeights(-1.0, 0.5, 0.5);
const Vec3 yWeights(0.0, -1.0, 1.0);
const Vec3 localPosition(0.4, 0.3, 0.2);
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new LocalCoordinatesSite(0, 1, 2, originWeights, xWeights, yWeights, localPosition));
CustomExternalForce* forceField = new CustomExternalForce("2*x^2+3*y^2+4*z^2");
system.addForce(forceField);
vector<double> params;
forceField->addParticle(0, params);
forceField->addParticle(1, params);
forceField->addParticle(2, params);
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4), positions2(4), positions3(4);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < 100; i++) {
// Set the particles at random positions.
Vec3 xdir, ydir, zdir;
do {
for (int j = 0; j < 3; j++)
positions[j] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
xdir = positions[0]*xWeights[0] + positions[1]*xWeights[1] + positions[2]*xWeights[2];
ydir = positions[0]*yWeights[0] + positions[1]*yWeights[1] + positions[2]*yWeights[2];
zdir = xdir.cross(ydir);
if (sqrt(xdir.dot(xdir)) > 0.1 && sqrt(ydir.dot(ydir)) > 0.1 && sqrt(zdir.dot(zdir)) > 0.1)
break; // These positions give a reasonable coordinate system.
} while (true);
context.setPositions(positions);
context.applyConstraints(0.0001);
// See if the virtual site is positioned correctly.
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 origin = pos[0]*originWeights[0] + pos[1]*originWeights[1] + pos[2]*originWeights[2];
xdir /= sqrt(xdir.dot(xdir));
zdir /= sqrt(zdir.dot(zdir));
ydir = zdir.cross(xdir);
ASSERT_EQUAL_VEC(origin+xdir*localPosition[0]+ydir*localPosition[1]+zdir*localPosition[2], pos[3], 1e-5);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
double norm = 0.0;
for (int i = 0; i < 3; ++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-2;
double step = 0.5*delta/norm;
for (int i = 0; i < 3; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
}
context.setPositions(positions2);
context.applyConstraints(0.0001);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
context.applyConstraints(0.0001);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/delta, 1e-3)
}
}
/**
* Make sure that energy, linear momentum, and angular momentum are all conserved
* when using virtual sites.
*/
void testConservationLaws() {
System system;
NonbondedForce* forceField = new NonbondedForce();
system.addForce(forceField);
vector<Vec3> positions;
// Create a linear molecule with a TwoParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.4, 0.6));
system.addConstraint(0, 1, 2.0);
for (int i = 0; i < 3; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i, j, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(2, 0, 0));
positions.push_back(Vec3());
// Create a planar molecule with a ThreeParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(6, new ThreeParticleAverageSite(3, 4, 5, 0.3, 0.5, 0.2));
system.addConstraint(3, 4, 1.0);
system.addConstraint(3, 5, 1.0);
system.addConstraint(4, 5, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+3, j+3, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 1));
positions.push_back(Vec3(1, 0, 1));
positions.push_back(Vec3(0, 1, 1));
positions.push_back(Vec3());
// Create a tetrahedral molecule with an OutOfPlane virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(10, new OutOfPlaneSite(7, 8, 9, 0.3, 0.5, 0.2));
system.addConstraint(7, 8, 1.0);
system.addConstraint(7, 9, 1.0);
system.addConstraint(8, 9, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+7, j+7, 0, 1, 0);
}
positions.push_back(Vec3(1, 0, -1));
positions.push_back(Vec3(2, 0, -1));
positions.push_back(Vec3(1, 1, -1));
positions.push_back(Vec3());
// Create a molecule with a LocalCoordinatesSite virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(14, new LocalCoordinatesSite(11, 12, 13, Vec3(0.3, 0.3, 0.4), Vec3(1.0, -0.5, -0.5), Vec3(0, -1.0, 1.0), Vec3(0.2, 0.2, 1.0)));
system.addConstraint(11, 12, 1.0);
system.addConstraint(11, 13, 1.0);
system.addConstraint(12, 13, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+11, j+11, 0, 1, 0);
}
positions.push_back(Vec3(1, 2, 0));
positions.push_back(Vec3(2, 2, 0));
positions.push_back(Vec3(1, 3, 0));
positions.push_back(Vec3());
// Simulate it and check conservation laws.
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
context.setPositions(positions);
context.applyConstraints(0.0001);
int numParticles = system.getNumParticles();
double initialEnergy;
Vec3 initialMomentum, initialAngularMomentum;
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
const vector<Vec3>& pos = state.getPositions();
const vector<Vec3>& vel = state.getVelocities();
const vector<Vec3>& f = state.getForces();
double energy = state.getPotentialEnergy();
for (int j = 0; j < numParticles; j++) {
Vec3 v = vel[j] + f[j]*0.5*integrator.getStepSize();
energy += 0.5*system.getParticleMass(j)*v.dot(v);
}
if (i == 0)
initialEnergy = energy;
else
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
Vec3 momentum;
for (int j = 0; j < numParticles; j++)
momentum += vel[j]*system.getParticleMass(j);
if (i == 0)
initialMomentum = momentum;
else
ASSERT_EQUAL_VEC(initialMomentum, momentum, 0.02);
Vec3 angularMomentum;
for (int j = 0; j < numParticles; j++)
angularMomentum += pos[j].cross(vel[j])*system.getParticleMass(j);
if (i == 0)
initialAngularMomentum = angularMomentum;
else
ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 0.03);
integrator.step(1);
}
}
#include "CudaTests.h"
#include "TestVirtualSites.h"
/**
* Make sure that atom reordering respects virtual sites.
......@@ -524,64 +138,6 @@ void testReordering() {
}
}
/**
* Test a System where multiple virtual sites are all calculated from the same particles.
*/
void testOverlappingSites() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
NonbondedForce* nonbonded = new NonbondedForce();
system.addForce(nonbonded);
nonbonded->addParticle(1.0, 0.0, 0.0);
nonbonded->addParticle(-0.5, 0.0, 0.0);
nonbonded->addParticle(-0.5, 0.0, 0.0);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(10, 0, 0));
positions.push_back(Vec3(0, 10, 0));
for (int i = 0; i < 20; i++) {
system.addParticle(0.0);
double u = 0.1*((i+1)%4);
double v = 0.05*i;
system.setVirtualSite(3+i, new ThreeParticleAverageSite(0, 1, 2, u, v, 1-u-v));
nonbonded->addParticle(i%2 == 0 ? -1.0 : 1.0, 0.0, 0.0);
positions.push_back(Vec3());
}
VerletIntegrator i1(0.002);
VerletIntegrator i2(0.002);
Context c1(system, i1, Platform::getPlatformByName("Reference"));
Context c2(system, i2, platform);
c1.setPositions(positions);
c2.setPositions(positions);
c1.applyConstraints(0.0001);
c2.applyConstraints(0.0001);
State s1 = c1.getState(State::Positions | State::Forces);
State s2 = c2.getState(State::Positions | State::Forces);
for (int i = 0; i < system.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getPositions()[i], s2.getPositions()[i], 1e-5);
for (int i = 0; i < 3; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testMasslessParticle();
testTwoParticleAverage();
testThreeParticleAverage();
testOutOfPlane();
testLocalCoordinates();
testConservationLaws();
testReordering();
testOverlappingSites();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
testReordering();
}
platforms/opencl/tests/TestOpenCLVariableLangevinIntegrator.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,311 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the OpenCL implementation of VariableLangevinIntegrator.
*/
#include "OpenCLTests.h"
#include "TestVariableLangevinIntegrator.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenCLPlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VariableLangevinIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
static OpenCLPlatform platform;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VariableLangevinIntegrator integrator(0, 0.1, 1e-6);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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 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 time = state.getTime();
double expectedDist = 1.5+0.5*std::exp(-0.05*time)*std::cos(freq*time);
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);
}
// Now 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);
double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
for (int i = 0; i < 1000; ++i) {
state = context.getState(State::Energy);
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
integrator.step(1);
}
}
void testTemperature() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 5.0, 5e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; ++i)
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Let it equilibrate.
integrator.step(5000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 5000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(5);
}
ke /= 5000;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.1);
}
void testConstraints() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 2.0, 1e-5);
integrator.setConstraintTolerance(1e-5);
integrator.setRandomNumberSeed(0);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numParticles-1; ++i)
system.addConstraint(i, i+1, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-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 < numParticles-1; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
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]));
ASSERT_EQUAL_TOL(1.0, dist, 2e-5);
}
integrator.step(1);
}
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VariableLangevinIntegrator integrator(300.0, 2.0, 0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
void testRandomSeed() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 2.0, 1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(0, 0, 0);
}
// Try twice with the same random seed.
integrator.setRandomNumberSeed(5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state1 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state2 = context.getState(State::Positions);
// Try twice with a different random seed.
integrator.setRandomNumberSeed(10);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state3 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state4 = context.getState(State::Positions);
// Compare the results.
for (int i = 0; i < numParticles; i++) {
for (int j = 0; j < 3; j++) {
ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
}
}
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
int numParticles = 0;
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
++numParticles;
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableLangevinIntegrator integrator(temp, 6.0, 1e-4);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(2.0);
// Make sure the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 1000; ++i) {
double t = 2.0 + 0.02 * (i + 1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
}
ke /= 1000;
double expected = 1.5 * numParticles * BOLTZ * temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.01);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
testSingleBond();
testTemperature();
testConstraints();
testConstrainedMasslessParticles();
testRandomSeed();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/opencl/tests/TestOpenCLVariableVerletIntegrator.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,289 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the OpenCL implementation of VariableVerletIntegrator.
*/
#include "OpenCLTests.h"
#include "TestVariableVerletIntegrator.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenCLPlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VariableVerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
static OpenCLPlatform platform;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VariableVerletIntegrator integrator(1e-6);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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();
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);
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);
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();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
integrator.step(1);
}
}
void testConstraints() {
const int numParticles = 8;
const int numConstraints = 5;
System system;
VariableVerletIntegrator integrator(1e-5);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(1, 2, 1.0);
system.addConstraint(2, 3, 1.0);
system.addConstraint(4, 5, 1.0);
system.addConstraint(6, 7, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
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]));
ASSERT_EQUAL_TOL(distance, dist, 1e-4);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
double finalTime = context.getState(State::Positions).getTime();
ASSERT(finalTime > 0.1);
// Now try the stepTo() method.
finalTime += 0.5;
integrator.stepTo(finalTime);
ASSERT_EQUAL(finalTime, context.getState(State::Positions).getTime());
}
void testConstrainedClusters() {
const int numParticles = 7;
System system;
VariableVerletIntegrator integrator(1e-5);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i > 1 ? 1.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(0, 2, 1.0);
system.addConstraint(0, 3, 1.0);
system.addConstraint(0, 4, 1.0);
system.addConstraint(1, 5, 1.0);
system.addConstraint(1, 6, 1.0);
system.addConstraint(2, 3, sqrt(2.0));
system.addConstraint(2, 4, sqrt(2.0));
system.addConstraint(3, 4, sqrt(2.0));
system.addConstraint(5, 6, sqrt(2.0));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(-1, 0, 0);
positions[3] = Vec3(0, 1, 0);
positions[4] = Vec3(0, 0, 1);
positions[5] = Vec3(2, 0, 0);
positions[6] = Vec3(1, 1, 0);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i)
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
for (int j = 0; j < system.getNumConstraints(); ++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]));
ASSERT_EQUAL_TOL(distance, dist, 2e-5);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
ASSERT(context.getState(State::Positions).getTime() > 0.1);
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VariableVerletIntegrator integrator(0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableVerletIntegrator integrator(1e-5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(1.0);
// Simulate it and see whether energy remains constant.
State state0 = context.getState(State::Energy);
double initialEnergy = state0.getKineticEnergy() + state0.getPotentialEnergy();
for (int i = 0; i < 20; i++) {
double t = 1.0 + 0.05*(i+1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
double energy = state.getKineticEnergy() + state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
}
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
testSingleBond();
testConstraints();
testConstrainedClusters();
testConstrainedMasslessParticles();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/opencl/tests/TestOpenCLVerletIntegrator.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,223 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the OpenCL implementation of VerletIntegrator.
*/
#include "OpenCLTests.h"
#include "TestVerletIntegrator.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenCLPlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
static OpenCLPlatform platform;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VerletIntegrator integrator(0.01);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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();
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);
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);
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();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
void testConstraints() {
const int numParticles = 8;
const int numConstraints = 5;
System system;
VerletIntegrator integrator(0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(1, 2, 1.0);
system.addConstraint(2, 3, 1.0);
system.addConstraint(4, 5, 1.0);
system.addConstraint(6, 7, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
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]));
ASSERT_EQUAL_TOL(distance, dist, 1e-4);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
void runPlatformTests() {
}
void testConstrainedClusters() {
const int numParticles = 7;
System system;
VerletIntegrator integrator(0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i > 1 ? 1.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(0, 2, 1.0);
system.addConstraint(0, 3, 1.0);
system.addConstraint(0, 4, 1.0);
system.addConstraint(1, 5, 1.0);
system.addConstraint(1, 6, 1.0);
system.addConstraint(2, 3, sqrt(2.0));
system.addConstraint(2, 4, sqrt(2.0));
system.addConstraint(3, 4, sqrt(2.0));
system.addConstraint(5, 6, sqrt(2.0));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(-1, 0, 0);
positions[3] = Vec3(0, 1, 0);
positions[4] = Vec3(0, 0, 1);
positions[5] = Vec3(2, 0, 0);
positions[6] = Vec3(1, 1, 0);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i)
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
for (int j = 0; j < system.getNumConstraints(); ++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]));
ASSERT_EQUAL_TOL(distance, dist, 2e-5);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VerletIntegrator integrator(0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
testSingleBond();
testConstraints();
testConstrainedClusters();
testConstrainedMasslessParticles();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/opencl/tests/TestOpenCLVirtualSites.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2012-2014 Stanford University and the Authors. *
* Portions copyright (c) 2012-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,394 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the OpenCL implementation of virtual sites.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenCLPlatform.h"
#include "openmm/CustomBondForce.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/VirtualSite.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
static OpenCLPlatform platform;
/**
* Check that massless particles are handled correctly.
*/
void testMasslessParticle() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
CustomBondForce* bonds = new CustomBondForce("-1/r");
system.addForce(bonds);
vector<double> params;
bonds->addBond(0, 1, params);
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
vector<Vec3> velocities(2);
velocities[0] = Vec3(0, 0, 0);
velocities[1] = Vec3(0, 1, 0);
context.setVelocities(velocities);
// The second particle should move in a circular orbit around the first one.
// Compare it to the analytical solution.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Velocities | State::Forces);
double time = state.getTime();
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getPositions()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getVelocities()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(cos(time), sin(time), 0), state.getPositions()[1], 0.01);
ASSERT_EQUAL_VEC(Vec3(-sin(time), cos(time), 0), state.getVelocities()[1], 0.01);
integrator.step(1);
}
}
/**
* Test a TwoParticleAverageSite virtual site.
*/
void testTwoParticleAverage() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.8, 0.2));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.8+pos[1]*0.2, pos[2], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.1+0.3*0.8, 0, 0), state.getForces()[0], 1e-4);
ASSERT_EQUAL_VEC(Vec3(0.2+0.3*0.2, 0, 0), state.getForces()[1], 1e-4);
integrator.step(1);
}
}
/**
* Test a ThreeParticleAverageSite virtual site.
*/
void testThreeParticleAverage() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new ThreeParticleAverageSite(0, 1, 2, 0.2, 0.3, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.2+pos[1]*0.3+pos[2]*0.5, pos[3], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4*0.2, 0, 0), state.getForces()[0], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.2+0.4*0.3, 0, 0), state.getForces()[1], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.3+0.4*0.5, 0, 0), state.getForces()[2], 1e-5);
integrator.step(1);
}
}
/**
* Test an OutOfPlaneSite virtual site.
*/
void testOutOfPlane() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new OutOfPlaneSite(0, 1, 2, 0.3, 0.4, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 v12 = pos[1]-pos[0];
Vec3 v13 = pos[2]-pos[0];
Vec3 cross = v12.cross(v13);
ASSERT_EQUAL_VEC(pos[0]+v12*0.3+v13*0.4+cross*0.5, pos[3], 1e-5);
const vector<Vec3>& f = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0.1+0.2+0.3+0.4, 0, 0), f[0]+f[1]+f[2], 1e-5);
Vec3 f2(0.4*0.3, 0.4*0.5*v13[2], -0.4*0.5*v13[1]);
Vec3 f3(0.4*0.4, -0.4*0.5*v12[2], 0.4*0.5*v12[1]);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4, 0, 0)-f2-f3, f[0], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.2, 0, 0)+f2, f[1], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.3, 0, 0)+f3, f[2], 1e-5);
integrator.step(1);
}
}
/**
* Test a LocalCoordinatesSite virtual site.
*/
void testLocalCoordinates() {
const Vec3 originWeights(0.2, 0.3, 0.5);
const Vec3 xWeights(-1.0, 0.5, 0.5);
const Vec3 yWeights(0.0, -1.0, 1.0);
const Vec3 localPosition(0.4, 0.3, 0.2);
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new LocalCoordinatesSite(0, 1, 2, originWeights, xWeights, yWeights, localPosition));
CustomExternalForce* forceField = new CustomExternalForce("2*x^2+3*y^2+4*z^2");
system.addForce(forceField);
vector<double> params;
forceField->addParticle(0, params);
forceField->addParticle(1, params);
forceField->addParticle(2, params);
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4), positions2(4), positions3(4);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < 100; i++) {
// Set the particles at random positions.
Vec3 xdir, ydir, zdir;
do {
for (int j = 0; j < 3; j++)
positions[j] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
xdir = positions[0]*xWeights[0] + positions[1]*xWeights[1] + positions[2]*xWeights[2];
ydir = positions[0]*yWeights[0] + positions[1]*yWeights[1] + positions[2]*yWeights[2];
zdir = xdir.cross(ydir);
if (sqrt(xdir.dot(xdir)) > 0.1 && sqrt(ydir.dot(ydir)) > 0.1 && sqrt(zdir.dot(zdir)) > 0.1)
break; // These positions give a reasonable coordinate system.
} while (true);
context.setPositions(positions);
context.applyConstraints(0.0001);
// See if the virtual site is positioned correctly.
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 origin = pos[0]*originWeights[0] + pos[1]*originWeights[1] + pos[2]*originWeights[2];
xdir /= sqrt(xdir.dot(xdir));
zdir /= sqrt(zdir.dot(zdir));
ydir = zdir.cross(xdir);
ASSERT_EQUAL_VEC(origin+xdir*localPosition[0]+ydir*localPosition[1]+zdir*localPosition[2], pos[3], 1e-5);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
double norm = 0.0;
for (int i = 0; i < 3; ++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-2;
double step = 0.5*delta/norm;
for (int i = 0; i < 3; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
}
context.setPositions(positions2);
context.applyConstraints(0.0001);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
context.applyConstraints(0.0001);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/delta, 1e-3)
}
}
/**
* Make sure that energy, linear momentum, and angular momentum are all conserved
* when using virtual sites.
*/
void testConservationLaws() {
System system;
NonbondedForce* forceField = new NonbondedForce();
system.addForce(forceField);
vector<Vec3> positions;
// Create a linear molecule with a TwoParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.4, 0.6));
system.addConstraint(0, 1, 2.0);
for (int i = 0; i < 3; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i, j, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(2, 0, 0));
positions.push_back(Vec3());
// Create a planar molecule with a ThreeParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(6, new ThreeParticleAverageSite(3, 4, 5, 0.3, 0.5, 0.2));
system.addConstraint(3, 4, 1.0);
system.addConstraint(3, 5, 1.0);
system.addConstraint(4, 5, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+3, j+3, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 1));
positions.push_back(Vec3(1, 0, 1));
positions.push_back(Vec3(0, 1, 1));
positions.push_back(Vec3());
// Create a tetrahedral molecule with an OutOfPlane virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(10, new OutOfPlaneSite(7, 8, 9, 0.3, 0.5, 0.2));
system.addConstraint(7, 8, 1.0);
system.addConstraint(7, 9, 1.0);
system.addConstraint(8, 9, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+7, j+7, 0, 1, 0);
}
positions.push_back(Vec3(1, 0, -1));
positions.push_back(Vec3(2, 0, -1));
positions.push_back(Vec3(1, 1, -1));
positions.push_back(Vec3());
// Create a molecule with a LocalCoordinatesSite virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(14, new LocalCoordinatesSite(11, 12, 13, Vec3(0.3, 0.3, 0.4), Vec3(1.0, -0.5, -0.5), Vec3(0, -1.0, 1.0), Vec3(0.2, 0.2, 1.0)));
system.addConstraint(11, 12, 1.0);
system.addConstraint(11, 13, 1.0);
system.addConstraint(12, 13, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+11, j+11, 0, 1, 0);
}
positions.push_back(Vec3(1, 2, 0));
positions.push_back(Vec3(2, 2, 0));
positions.push_back(Vec3(1, 3, 0));
positions.push_back(Vec3());
// Simulate it and check conservation laws.
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
context.setPositions(positions);
context.applyConstraints(0.0001);
int numParticles = system.getNumParticles();
double initialEnergy;
Vec3 initialMomentum, initialAngularMomentum;
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
const vector<Vec3>& pos = state.getPositions();
const vector<Vec3>& vel = state.getVelocities();
const vector<Vec3>& f = state.getForces();
double energy = state.getPotentialEnergy();
for (int j = 0; j < numParticles; j++) {
Vec3 v = vel[j] + f[j]*0.5*integrator.getStepSize();
energy += 0.5*system.getParticleMass(j)*v.dot(v);
}
if (i == 0)
initialEnergy = energy;
else
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
Vec3 momentum;
for (int j = 0; j < numParticles; j++)
momentum += vel[j]*system.getParticleMass(j);
if (i == 0)
initialMomentum = momentum;
else
ASSERT_EQUAL_VEC(initialMomentum, momentum, 0.02);
Vec3 angularMomentum;
for (int j = 0; j < numParticles; j++)
angularMomentum += pos[j].cross(vel[j])*system.getParticleMass(j);
if (i == 0)
initialAngularMomentum = angularMomentum;
else
ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 0.05);
integrator.step(1);
}
}
#include "OpenCLTests.h"
#include "TestVirtualSites.h"
/**
* Make sure that atom reordering respects virtual sites.
......@@ -524,64 +138,6 @@ void testReordering() {
}
}
/**
* Test a System where multiple virtual sites are all calculated from the same particles.
*/
void testOverlappingSites() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
NonbondedForce* nonbonded = new NonbondedForce();
system.addForce(nonbonded);
nonbonded->addParticle(1.0, 0.0, 0.0);
nonbonded->addParticle(-0.5, 0.0, 0.0);
nonbonded->addParticle(-0.5, 0.0, 0.0);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(10, 0, 0));
positions.push_back(Vec3(0, 10, 0));
for (int i = 0; i < 20; i++) {
system.addParticle(0.0);
double u = 0.1*((i+1)%4);
double v = 0.05*i;
system.setVirtualSite(3+i, new ThreeParticleAverageSite(0, 1, 2, u, v, 1-u-v));
nonbonded->addParticle(i%2 == 0 ? -1.0 : 1.0, 0.0, 0.0);
positions.push_back(Vec3());
}
VerletIntegrator i1(0.002);
VerletIntegrator i2(0.002);
Context c1(system, i1, Platform::getPlatformByName("Reference"));
Context c2(system, i2, platform);
c1.setPositions(positions);
c2.setPositions(positions);
c1.applyConstraints(0.0001);
c2.applyConstraints(0.0001);
State s1 = c1.getState(State::Positions | State::Forces);
State s2 = c2.getState(State::Positions | State::Forces);
for (int i = 0; i < system.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getPositions()[i], s2.getPositions()[i], 1e-5);
for (int i = 0; i < 3; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
testMasslessParticle();
testTwoParticleAverage();
testThreeParticleAverage();
testOutOfPlane();
testLocalCoordinates();
testConservationLaws();
testReordering();
testOverlappingSites();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
testReordering();
}
platforms/reference/tests/TestReferenceVariableLangevinIntegrator.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,309 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of VariableLangevinIntegrator.
*/
#include "ReferenceTests.h"
#include "TestVariableLangevinIntegrator.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VariableLangevinIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VariableLangevinIntegrator integrator(0, 0.1, 1e-6);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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 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 time = state.getTime();
double expectedDist = 1.5+0.5*std::exp(-0.05*time)*std::cos(freq*time);
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);
}
// Now 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);
double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
for (int i = 0; i < 1000; ++i) {
state = context.getState(State::Energy);
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
integrator.step(1);
}
}
void testTemperature() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 5.0, 5e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; ++i)
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Let it equilibrate.
integrator.step(5000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 5000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(5);
}
ke /= 5000;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.1);
}
void testConstraints() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 2.0, 1e-5);
integrator.setConstraintTolerance(1e-5);
integrator.setRandomNumberSeed(0);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numParticles-1; ++i)
system.addConstraint(i, i+1, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-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 < numParticles-1; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
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]));
ASSERT_EQUAL_TOL(1.0, dist, 2e-5);
}
integrator.step(1);
}
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VariableLangevinIntegrator integrator(300.0, 2.0, 0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities | State::Positions);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
void testRandomSeed() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 2.0, 1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(0, 0, 0);
}
// Try twice with the same random seed.
integrator.setRandomNumberSeed(5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state1 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state2 = context.getState(State::Positions);
// Try twice with a different random seed.
integrator.setRandomNumberSeed(10);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state3 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state4 = context.getState(State::Positions);
// Compare the results.
for (int i = 0; i < numParticles; i++) {
for (int j = 0; j < 3; j++) {
ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
}
}
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
int numParticles = 0;
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
++numParticles;
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableLangevinIntegrator integrator(temp, 6.0, 1e-4);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(2.0);
// Make sure the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 1000; ++i) {
double t = 2.0 + 0.02 * (i + 1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
}
ke /= 1000;
double expected = 1.5 * numParticles * BOLTZ * temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.01);
}
int main() {
try {
testSingleBond();
testTemperature();
testConstraints();
testConstrainedMasslessParticles();
testRandomSeed();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/reference/tests/TestReferenceVariableVerletIntegrator.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,286 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of VariableVerletIntegrator.
*/
#include "ReferenceTests.h"
#include "TestVariableVerletIntegrator.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VariableVerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VariableVerletIntegrator integrator(1e-6);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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();
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);
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);
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();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
integrator.step(1);
}
ASSERT(state.getTime() > 1.0);
}
void testConstraints() {
const int numParticles = 8;
const double temp = 500.0;
System system;
VariableVerletIntegrator integrator(1e-5);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i%2 == 0 ? 5.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numParticles-1; ++i)
system.addConstraint(i, i+1, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
for (int j = 0; j < system.getNumConstraints(); ++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]));
ASSERT_EQUAL_TOL(distance, dist, 2e-5);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
double finalTime = context.getState(State::Positions).getTime();
ASSERT(finalTime > 0.1);
// Now try the stepTo() method.
finalTime += 0.5;
integrator.stepTo(finalTime);
ASSERT_EQUAL(finalTime, context.getState(State::Positions).getTime());
}
void testConstrainedClusters() {
const int numParticles = 7;
const double temp = 500.0;
System system;
VariableVerletIntegrator integrator(1e-5);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i > 1 ? 1.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(0, 2, 1.0);
system.addConstraint(0, 3, 1.0);
system.addConstraint(0, 4, 1.0);
system.addConstraint(1, 5, 1.0);
system.addConstraint(1, 6, 1.0);
system.addConstraint(2, 3, sqrt(2.0));
system.addConstraint(2, 4, sqrt(2.0));
system.addConstraint(3, 4, sqrt(2.0));
system.addConstraint(5, 6, sqrt(2.0));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(-1, 0, 0);
positions[3] = Vec3(0, 1, 0);
positions[4] = Vec3(0, 0, 1);
positions[5] = Vec3(2, 0, 0);
positions[6] = Vec3(1, 1, 0);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i)
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
for (int j = 0; j < system.getNumConstraints(); ++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]));
ASSERT_EQUAL_TOL(distance, dist, 2e-5);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
ASSERT(context.getState(State::Positions).getTime() > 0.1);
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VariableVerletIntegrator integrator(0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities | State::Positions);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
void runPlatformTests() {
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableVerletIntegrator integrator(1e-5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(1.0);
// Simulate it and see whether energy remains constant.
State state0 = context.getState(State::Energy);
double initialEnergy = state0.getKineticEnergy() + state0.getPotentialEnergy();
for (int i = 0; i < 20; i++) {
double t = 1.0 + 0.05*(i+1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
double energy = state.getKineticEnergy() + state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
}
}
int main() {
try {
testSingleBond();
testConstraints();
testConstrainedClusters();
testConstrainedMasslessParticles();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/reference/tests/TestReferenceVerletIntegrator.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,218 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of VerletIntegrator.
*/
#include "ReferenceTests.h"
#include "TestVerletIntegrator.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VerletIntegrator integrator(0.01);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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();
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);
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);
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();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
void testConstraints() {
const int numParticles = 8;
const double temp = 500.0;
System system;
VerletIntegrator integrator(0.002);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i%2 == 0 ? 5.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numParticles-1; ++i)
system.addConstraint(i, i+1, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
for (int j = 0; j < system.getNumConstraints(); ++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]));
ASSERT_EQUAL_TOL(distance, dist, 2e-5);
}
double energy = state.getPotentialEnergy()+state.getKineticEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
void testConstrainedClusters() {
const int numParticles = 7;
const double temp = 500.0;
System system;
VerletIntegrator integrator(0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i > 1 ? 1.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(0, 2, 1.0);
system.addConstraint(0, 3, 1.0);
system.addConstraint(0, 4, 1.0);
system.addConstraint(1, 5, 1.0);
system.addConstraint(1, 6, 1.0);
system.addConstraint(2, 3, sqrt(2.0));
system.addConstraint(2, 4, sqrt(2.0));
system.addConstraint(3, 4, sqrt(2.0));
system.addConstraint(5, 6, sqrt(2.0));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(-1, 0, 0);
positions[3] = Vec3(0, 1, 0);
positions[4] = Vec3(0, 0, 1);
positions[5] = Vec3(2, 0, 0);
positions[6] = Vec3(1, 1, 0);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i)
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
for (int j = 0; j < system.getNumConstraints(); ++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]));
ASSERT_EQUAL_TOL(distance, dist, 2e-5);
}
double energy = state.getPotentialEnergy()+state.getKineticEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VerletIntegrator integrator(0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities | State::Positions);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
int main() {
try {
testSingleBond();
testConstraints();
testConstrainedClusters();
testConstrainedMasslessParticles();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/reference/tests/TestReferenceVirtualSites.cpp
View file @ 8b5d3dc0
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2012-2014 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,408 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of virtual sites.
*/
#include "ReferenceTests.h"
#include "TestVirtualSites.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/CustomBondForce.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/VirtualSite.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
ReferencePlatform platform;
/**
* Check that massless particles are handled correctly.
*/
void testMasslessParticle() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
CustomBondForce* bonds = new CustomBondForce("-1/r");
system.addForce(bonds);
vector<double> params;
bonds->addBond(0, 1, params);
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
vector<Vec3> velocities(2);
velocities[0] = Vec3(0, 0, 0);
velocities[1] = Vec3(0, 1, 0);
context.setVelocities(velocities);
// The second particle should move in a circular orbit around the first one.
// Compare it to the analytical solution.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Velocities | State::Forces);
double time = state.getTime();
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getPositions()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getVelocities()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(cos(time), sin(time), 0), state.getPositions()[1], 0.01);
ASSERT_EQUAL_VEC(Vec3(-sin(time), cos(time), 0), state.getVelocities()[1], 0.01);
integrator.step(1);
}
}
/**
* Test a TwoParticleAverageSite virtual site.
*/
void testTwoParticleAverage() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.8, 0.2));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.8+pos[1]*0.2, pos[2], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.1+0.3*0.8, 0, 0), state.getForces()[0], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.2+0.3*0.2, 0, 0), state.getForces()[1], 1e-10);
integrator.step(1);
}
}
/**
* Test a ThreeParticleAverageSite virtual site.
*/
void testThreeParticleAverage() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new ThreeParticleAverageSite(0, 1, 2, 0.2, 0.3, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.2+pos[1]*0.3+pos[2]*0.5, pos[3], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4*0.2, 0, 0), state.getForces()[0], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.2+0.4*0.3, 0, 0), state.getForces()[1], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.3+0.4*0.5, 0, 0), state.getForces()[2], 1e-10);
integrator.step(1);
}
}
/**
* Test an OutOfPlaneSite virtual site.
*/
void testOutOfPlane() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new OutOfPlaneSite(0, 1, 2, 0.3, 0.4, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 v12 = pos[1]-pos[0];
Vec3 v13 = pos[2]-pos[0];
Vec3 cross = v12.cross(v13);
ASSERT_EQUAL_VEC(pos[0]+v12*0.3+v13*0.4+cross*0.5, pos[3], 1e-10);
const vector<Vec3>& f = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0.1+0.2+0.3+0.4, 0, 0), f[0]+f[1]+f[2], 1e-10);
Vec3 f2(0.4*0.3, 0.4*0.5*v13[2], -0.4*0.5*v13[1]);
Vec3 f3(0.4*0.4, -0.4*0.5*v12[2], 0.4*0.5*v12[1]);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4, 0, 0)-f2-f3, f[0], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.2, 0, 0)+f2, f[1], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.3, 0, 0)+f3, f[2], 1e-10);
integrator.step(1);
}
}
/**
* Test a LocalCoordinatesSite virtual site.
*/
void testLocalCoordinates() {
const Vec3 originWeights(0.2, 0.3, 0.5);
const Vec3 xWeights(-1.0, 0.5, 0.5);
const Vec3 yWeights(0.0, -1.0, 1.0);
const Vec3 localPosition(0.4, 0.3, 0.2);
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new LocalCoordinatesSite(0, 1, 2, originWeights, xWeights, yWeights, localPosition));
CustomExternalForce* forceField = new CustomExternalForce("2*x^2+3*y^2+4*z^2");
system.addForce(forceField);
vector<double> params;
forceField->addParticle(0, params);
forceField->addParticle(1, params);
forceField->addParticle(2, params);
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4), positions2(4), positions3(4);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < 100; i++) {
// Set the particles at random positions.
Vec3 xdir, ydir, zdir;
do {
for (int j = 0; j < 3; j++)
positions[j] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
xdir = positions[0]*xWeights[0] + positions[1]*xWeights[1] + positions[2]*xWeights[2];
ydir = positions[0]*yWeights[0] + positions[1]*yWeights[1] + positions[2]*yWeights[2];
zdir = xdir.cross(ydir);
if (sqrt(xdir.dot(xdir)) > 0.1 && sqrt(ydir.dot(ydir)) > 0.1 && sqrt(zdir.dot(zdir)) > 0.1)
break; // These positions give a reasonable coordinate system.
} while (true);
context.setPositions(positions);
context.applyConstraints(0.0001);
// See if the virtual site is positioned correctly.
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 origin = pos[0]*originWeights[0] + pos[1]*originWeights[1] + pos[2]*originWeights[2];
xdir /= sqrt(xdir.dot(xdir));
zdir /= sqrt(zdir.dot(zdir));
ydir = zdir.cross(xdir);
ASSERT_EQUAL_VEC(origin+xdir*localPosition[0]+ydir*localPosition[1]+zdir*localPosition[2], pos[3], 1e-10);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
double norm = 0.0;
for (int i = 0; i < 3; ++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-2;
double step = 0.5*delta/norm;
for (int i = 0; i < 3; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
}
context.setPositions(positions2);
context.applyConstraints(0.0001);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
context.applyConstraints(0.0001);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/delta, 1e-3)
}
}
/**
* Make sure that energy, linear momentum, and angular momentum are all conserved
* when using virtual sites.
*/
void testConservationLaws() {
System system;
NonbondedForce* forceField = new NonbondedForce();
system.addForce(forceField);
vector<Vec3> positions;
// Create a linear molecule with a TwoParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.4, 0.6));
system.addConstraint(0, 1, 2.0);
for (int i = 0; i < 3; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i, j, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(2, 0, 0));
positions.push_back(Vec3());
// Create a planar molecule with a ThreeParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(6, new ThreeParticleAverageSite(3, 4, 5, 0.3, 0.5, 0.2));
system.addConstraint(3, 4, 1.0);
system.addConstraint(3, 5, 1.0);
system.addConstraint(4, 5, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+3, j+3, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 1));
positions.push_back(Vec3(1, 0, 1));
positions.push_back(Vec3(0, 1, 1));
positions.push_back(Vec3());
// Create a tetrahedral molecule with an OutOfPlane virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(10, new OutOfPlaneSite(7, 8, 9, 0.3, 0.5, 0.2));
system.addConstraint(7, 8, 1.0);
system.addConstraint(7, 9, 1.0);
system.addConstraint(8, 9, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+7, j+7, 0, 1, 0);
}
positions.push_back(Vec3(1, 0, -1));
positions.push_back(Vec3(2, 0, -1));
positions.push_back(Vec3(1, 1, -1));
positions.push_back(Vec3());
// Create a molecule with a LocalCoordinatesSite virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(14, new LocalCoordinatesSite(11, 12, 13, Vec3(0.3, 0.3, 0.4), Vec3(1.0, -0.5, -0.5), Vec3(0, -1.0, 1.0), Vec3(0.2, 0.2, 1.0)));
system.addConstraint(11, 12, 1.0);
system.addConstraint(11, 13, 1.0);
system.addConstraint(12, 13, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+11, j+11, 0, 1, 0);
}
positions.push_back(Vec3(1, 2, 0));
positions.push_back(Vec3(2, 2, 0));
positions.push_back(Vec3(1, 3, 0));
positions.push_back(Vec3());
// Simulate it and check conservation laws.
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
context.setPositions(positions);
context.applyConstraints(0.0001);
int numParticles = system.getNumParticles();
double initialEnergy;
Vec3 initialMomentum, initialAngularMomentum;
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
const vector<Vec3>& pos = state.getPositions();
const vector<Vec3>& vel = state.getVelocities();
const vector<Vec3>& f = state.getForces();
double energy = state.getPotentialEnergy();
for (int j = 0; j < numParticles; j++) {
Vec3 v = vel[j] + f[j]*0.5*integrator.getStepSize();
energy += 0.5*system.getParticleMass(j)*v.dot(v);
}
if (i == 0)
initialEnergy = energy;
else
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
Vec3 momentum;
for (int j = 0; j < numParticles; j++)
momentum += vel[j]*system.getParticleMass(j);
if (i == 0)
initialMomentum = momentum;
else
ASSERT_EQUAL_VEC(initialMomentum, momentum, 1e-10);
Vec3 angularMomentum;
for (int j = 0; j < numParticles; j++)
angularMomentum += pos[j].cross(vel[j])*system.getParticleMass(j);
if (i == 0)
initialAngularMomentum = angularMomentum;
else
ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 1e-10);
integrator.step(1);
}
}
int main() {
try {
testMasslessParticle();
testTwoParticleAverage();
testThreeParticleAverage();
testOutOfPlane();
testLocalCoordinates();
testConservationLaws();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
tests/TestVariableLangevinIntegrator.h 0 → 100644
View file @ 8b5d3dc0
/* -------------------------------------------------------------------------- *
* 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) 2008-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VariableLangevinIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VariableLangevinIntegrator integrator(0, 0.1, 1e-6);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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 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 time = state.getTime();
double expectedDist = 1.5+0.5*std::exp(-0.05*time)*std::cos(freq*time);
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);
}
// Now 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);
double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
for (int i = 0; i < 1000; ++i) {
state = context.getState(State::Energy);
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
integrator.step(1);
}
}
void testTemperature() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 5.0, 5e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; ++i)
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Let it equilibrate.
integrator.step(5000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 5000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(5);
}
ke /= 5000;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.1);
}
void testConstraints() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 2.0, 1e-5);
integrator.setConstraintTolerance(1e-5);
integrator.setRandomNumberSeed(0);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numParticles-1; ++i)
system.addConstraint(i, i+1, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-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 < numParticles-1; ++j) {
Vec3 p1 = state.getPositions()[j];
Vec3 p2 = state.getPositions()[j+1];
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]));
ASSERT_EQUAL_TOL(1.0, dist, 2e-5);
}
integrator.step(1);
}
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VariableLangevinIntegrator integrator(300.0, 2.0, 0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
void testRandomSeed() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 2.0, 1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
}
system.addForce(forceField);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
velocities[i] = Vec3(0, 0, 0);
}
// Try twice with the same random seed.
integrator.setRandomNumberSeed(5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state1 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state2 = context.getState(State::Positions);
// Try twice with a different random seed.
integrator.setRandomNumberSeed(10);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state3 = context.getState(State::Positions);
context.reinitialize();
context.setPositions(positions);
context.setVelocities(velocities);
integrator.step(10);
State state4 = context.getState(State::Positions);
// Compare the results.
for (int i = 0; i < numParticles; i++) {
for (int j = 0; j < 3; j++) {
ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
}
}
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
int numParticles = 0;
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
++numParticles;
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableLangevinIntegrator integrator(temp, 6.0, 1e-4);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(2.0);
// Make sure the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 1000; ++i) {
double t = 2.0 + 0.02 * (i + 1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
}
ke /= 1000;
double expected = 1.5 * numParticles * BOLTZ * temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.01);
}
void runPlatformTests();
int main(int argc, char* argv[]) {
try {
initializeTests(argc, argv);
testSingleBond();
testTemperature();
testConstraints();
testConstrainedMasslessParticles();
testRandomSeed();
testArgonBox();
runPlatformTests();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
tests/TestVariableVerletIntegrator.h 0 → 100644
View file @ 8b5d3dc0
/* -------------------------------------------------------------------------- *
* 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) 2008-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VariableVerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VariableVerletIntegrator integrator(1e-6);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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();
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);
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);
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();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
integrator.step(1);
}
ASSERT(state.getTime() > 1.0);
}
void testConstraints() {
const int numParticles = 8;
const int numConstraints = 5;
System system;
VariableVerletIntegrator integrator(1e-5);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i%2 == 0 ? 5.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(1, 2, 1.0);
system.addConstraint(2, 3, 1.0);
system.addConstraint(4, 5, 1.0);
system.addConstraint(6, 7, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
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]));
ASSERT_EQUAL_TOL(distance, dist, 1e-4);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
double finalTime = context.getState(State::Positions).getTime();
ASSERT(finalTime > 0.1);
// Now try the stepTo() method.
finalTime += 0.5;
integrator.stepTo(finalTime);
ASSERT_EQUAL(finalTime, context.getState(State::Positions).getTime());
}
void testConstrainedClusters() {
const int numParticles = 7;
System system;
VariableVerletIntegrator integrator(1e-5);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i > 1 ? 1.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(0, 2, 1.0);
system.addConstraint(0, 3, 1.0);
system.addConstraint(0, 4, 1.0);
system.addConstraint(1, 5, 1.0);
system.addConstraint(1, 6, 1.0);
system.addConstraint(2, 3, sqrt(2.0));
system.addConstraint(2, 4, sqrt(2.0));
system.addConstraint(3, 4, sqrt(2.0));
system.addConstraint(5, 6, sqrt(2.0));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(-1, 0, 0);
positions[3] = Vec3(0, 1, 0);
positions[4] = Vec3(0, 0, 1);
positions[5] = Vec3(2, 0, 0);
positions[6] = Vec3(1, 1, 0);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i)
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
for (int j = 0; j < system.getNumConstraints(); ++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]));
ASSERT_EQUAL_TOL(distance, dist, 2e-5);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
ASSERT(context.getState(State::Positions).getTime() > 0.1);
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VariableVerletIntegrator integrator(0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableVerletIntegrator integrator(1e-5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(1.0);
// Simulate it and see whether energy remains constant.
State state0 = context.getState(State::Energy);
double initialEnergy = state0.getKineticEnergy() + state0.getPotentialEnergy();
for (int i = 0; i < 20; i++) {
double t = 1.0 + 0.05*(i+1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
double energy = state.getKineticEnergy() + state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
}
}
void runPlatformTests();
int main(int argc, char* argv[]) {
try {
initializeTests(argc, argv);
testSingleBond();
testConstraints();
testConstrainedClusters();
testConstrainedMasslessParticles();
testArgonBox();
runPlatformTests();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
tests/TestVerletIntegrator.h 0 → 100644
View file @ 8b5d3dc0
/* -------------------------------------------------------------------------- *
* 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) 2008-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VerletIntegrator integrator(0.01);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
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();
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);
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);
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();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
ASSERT_EQUAL_TOL(10.0, context.getState(0).getTime(), 1e-5);
}
void testConstraints() {
const int numParticles = 8;
const int numConstraints = 5;
System system;
VerletIntegrator integrator(0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i%2 == 0 ? 5.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(1, 2, 1.0);
system.addConstraint(2, 3, 1.0);
system.addConstraint(4, 5, 1.0);
system.addConstraint(6, 7, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
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]));
ASSERT_EQUAL_TOL(distance, dist, 1e-4);
}
double energy = state.getPotentialEnergy()+state.getKineticEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
void testConstrainedClusters() {
const int numParticles = 7;
System system;
VerletIntegrator integrator(0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i > 1 ? 1.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
system.addConstraint(0, 1, 1.0);
system.addConstraint(0, 2, 1.0);
system.addConstraint(0, 3, 1.0);
system.addConstraint(0, 4, 1.0);
system.addConstraint(1, 5, 1.0);
system.addConstraint(1, 6, 1.0);
system.addConstraint(2, 3, sqrt(2.0));
system.addConstraint(2, 4, sqrt(2.0));
system.addConstraint(3, 4, sqrt(2.0));
system.addConstraint(5, 6, sqrt(2.0));
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(-1, 0, 0);
positions[3] = Vec3(0, 1, 0);
positions[4] = Vec3(0, 0, 1);
positions[5] = Vec3(2, 0, 0);
positions[6] = Vec3(1, 1, 0);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i)
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy | State::Velocities | State::Forces);
for (int j = 0; j < system.getNumConstraints(); ++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]));
ASSERT_EQUAL_TOL(distance, dist, 2e-5);
}
double energy = state.getPotentialEnergy()+state.getKineticEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
void testConstrainedMasslessParticles() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
system.addConstraint(0, 1, 1.5);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VerletIntegrator integrator(0.01);
bool failed = false;
try {
// This should throw an exception.
Context context(system, integrator, platform);
}
catch (exception& ex) {
failed = true;
}
ASSERT(failed);
// Now make both particles massless, which should work.
system.setParticleMass(1, 0.0);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(300.0);
integrator.step(1);
State state = context.getState(State::Velocities);
ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
}
void runPlatformTests();
int main(int argc, char* argv[]) {
try {
initializeTests(argc, argv);
testSingleBond();
testConstraints();
testConstrainedClusters();
testConstrainedMasslessParticles();
runPlatformTests();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
tests/TestVirtualSites.h 0 → 100644
View file @ 8b5d3dc0
/* -------------------------------------------------------------------------- *
* 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) 2012-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/CustomBondForce.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/VirtualSite.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
/**
* Check that massless particles are handled correctly.
*/
void testMasslessParticle() {
System system;
system.addParticle(0.0);
system.addParticle(1.0);
CustomBondForce* bonds = new CustomBondForce("-1/r");
system.addForce(bonds);
vector<double> params;
bonds->addBond(0, 1, params);
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
vector<Vec3> velocities(2);
velocities[0] = Vec3(0, 0, 0);
velocities[1] = Vec3(0, 1, 0);
context.setVelocities(velocities);
// The second particle should move in a circular orbit around the first one.
// Compare it to the analytical solution.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Velocities | State::Forces);
double time = state.getTime();
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getPositions()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getVelocities()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(cos(time), sin(time), 0), state.getPositions()[1], 0.01);
ASSERT_EQUAL_VEC(Vec3(-sin(time), cos(time), 0), state.getVelocities()[1], 0.01);
integrator.step(1);
}
}
/**
* Test a TwoParticleAverageSite virtual site.
*/
void testTwoParticleAverage() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.8, 0.2));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.8+pos[1]*0.2, pos[2], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.1+0.3*0.8, 0, 0), state.getForces()[0], 1e-4);
ASSERT_EQUAL_VEC(Vec3(0.2+0.3*0.2, 0, 0), state.getForces()[1], 1e-4);
integrator.step(1);
}
}
/**
* Test a ThreeParticleAverageSite virtual site.
*/
void testThreeParticleAverage() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new ThreeParticleAverageSite(0, 1, 2, 0.2, 0.3, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.2+pos[1]*0.3+pos[2]*0.5, pos[3], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4*0.2, 0, 0), state.getForces()[0], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.2+0.4*0.3, 0, 0), state.getForces()[1], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.3+0.4*0.5, 0, 0), state.getForces()[2], 1e-5);
integrator.step(1);
}
}
/**
* Test an OutOfPlaneSite virtual site.
*/
void testOutOfPlane() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new OutOfPlaneSite(0, 1, 2, 0.3, 0.4, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 v12 = pos[1]-pos[0];
Vec3 v13 = pos[2]-pos[0];
Vec3 cross = v12.cross(v13);
ASSERT_EQUAL_VEC(pos[0]+v12*0.3+v13*0.4+cross*0.5, pos[3], 1e-5);
const vector<Vec3>& f = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0.1+0.2+0.3+0.4, 0, 0), f[0]+f[1]+f[2], 1e-5);
Vec3 f2(0.4*0.3, 0.4*0.5*v13[2], -0.4*0.5*v13[1]);
Vec3 f3(0.4*0.4, -0.4*0.5*v12[2], 0.4*0.5*v12[1]);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4, 0, 0)-f2-f3, f[0], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.2, 0, 0)+f2, f[1], 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.3, 0, 0)+f3, f[2], 1e-5);
integrator.step(1);
}
}
/**
* Test a LocalCoordinatesSite virtual site.
*/
void testLocalCoordinates() {
const Vec3 originWeights(0.2, 0.3, 0.5);
const Vec3 xWeights(-1.0, 0.5, 0.5);
const Vec3 yWeights(0.0, -1.0, 1.0);
const Vec3 localPosition(0.4, 0.3, 0.2);
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new LocalCoordinatesSite(0, 1, 2, originWeights, xWeights, yWeights, localPosition));
CustomExternalForce* forceField = new CustomExternalForce("2*x^2+3*y^2+4*z^2");
system.addForce(forceField);
vector<double> params;
forceField->addParticle(0, params);
forceField->addParticle(1, params);
forceField->addParticle(2, params);
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4), positions2(4), positions3(4);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < 100; i++) {
// Set the particles at random positions.
Vec3 xdir, ydir, zdir;
do {
for (int j = 0; j < 3; j++)
positions[j] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
xdir = positions[0]*xWeights[0] + positions[1]*xWeights[1] + positions[2]*xWeights[2];
ydir = positions[0]*yWeights[0] + positions[1]*yWeights[1] + positions[2]*yWeights[2];
zdir = xdir.cross(ydir);
if (sqrt(xdir.dot(xdir)) > 0.1 && sqrt(ydir.dot(ydir)) > 0.1 && sqrt(zdir.dot(zdir)) > 0.1)
break; // These positions give a reasonable coordinate system.
} while (true);
context.setPositions(positions);
context.applyConstraints(0.0001);
// See if the virtual site is positioned correctly.
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 origin = pos[0]*originWeights[0] + pos[1]*originWeights[1] + pos[2]*originWeights[2];
xdir /= sqrt(xdir.dot(xdir));
zdir /= sqrt(zdir.dot(zdir));
ydir = zdir.cross(xdir);
ASSERT_EQUAL_VEC(origin+xdir*localPosition[0]+ydir*localPosition[1]+zdir*localPosition[2], pos[3], 1e-5);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
double norm = 0.0;
for (int i = 0; i < 3; ++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-2;
double step = 0.5*delta/norm;
for (int i = 0; i < 3; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
}
context.setPositions(positions2);
context.applyConstraints(0.0001);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
context.applyConstraints(0.0001);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/delta, 1e-3)
}
}
/**
* Make sure that energy, linear momentum, and angular momentum are all conserved
* when using virtual sites.
*/
void testConservationLaws() {
System system;
NonbondedForce* forceField = new NonbondedForce();
system.addForce(forceField);
vector<Vec3> positions;
// Create a linear molecule with a TwoParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.4, 0.6));
system.addConstraint(0, 1, 2.0);
for (int i = 0; i < 3; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i, j, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(2, 0, 0));
positions.push_back(Vec3());
// Create a planar molecule with a ThreeParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(6, new ThreeParticleAverageSite(3, 4, 5, 0.3, 0.5, 0.2));
system.addConstraint(3, 4, 1.0);
system.addConstraint(3, 5, 1.0);
system.addConstraint(4, 5, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+3, j+3, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 1));
positions.push_back(Vec3(1, 0, 1));
positions.push_back(Vec3(0, 1, 1));
positions.push_back(Vec3());
// Create a tetrahedral molecule with an OutOfPlane virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(10, new OutOfPlaneSite(7, 8, 9, 0.3, 0.5, 0.2));
system.addConstraint(7, 8, 1.0);
system.addConstraint(7, 9, 1.0);
system.addConstraint(8, 9, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+7, j+7, 0, 1, 0);
}
positions.push_back(Vec3(1, 0, -1));
positions.push_back(Vec3(2, 0, -1));
positions.push_back(Vec3(1, 1, -1));
positions.push_back(Vec3());
// Create a molecule with a LocalCoordinatesSite virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(14, new LocalCoordinatesSite(11, 12, 13, Vec3(0.3, 0.3, 0.4), Vec3(1.0, -0.5, -0.5), Vec3(0, -1.0, 1.0), Vec3(0.2, 0.2, 1.0)));
system.addConstraint(11, 12, 1.0);
system.addConstraint(11, 13, 1.0);
system.addConstraint(12, 13, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+11, j+11, 0, 1, 0);
}
positions.push_back(Vec3(1, 2, 0));
positions.push_back(Vec3(2, 2, 0));
positions.push_back(Vec3(1, 3, 0));
positions.push_back(Vec3());
// Simulate it and check conservation laws.
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
context.setPositions(positions);
context.applyConstraints(0.0001);
int numParticles = system.getNumParticles();
double initialEnergy;
Vec3 initialMomentum, initialAngularMomentum;
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
const vector<Vec3>& pos = state.getPositions();
const vector<Vec3>& vel = state.getVelocities();
const vector<Vec3>& f = state.getForces();
double energy = state.getPotentialEnergy();
for (int j = 0; j < numParticles; j++) {
Vec3 v = vel[j] + f[j]*0.5*integrator.getStepSize();
energy += 0.5*system.getParticleMass(j)*v.dot(v);
}
if (i == 0)
initialEnergy = energy;
else
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
Vec3 momentum;
for (int j = 0; j < numParticles; j++)
momentum += vel[j]*system.getParticleMass(j);
if (i == 0)
initialMomentum = momentum;
else
ASSERT_EQUAL_VEC(initialMomentum, momentum, 0.02);
Vec3 angularMomentum;
for (int j = 0; j < numParticles; j++)
angularMomentum += pos[j].cross(vel[j])*system.getParticleMass(j);
if (i == 0)
initialAngularMomentum = angularMomentum;
else
ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 0.05);
integrator.step(1);
}
}
/**
* Test a System where multiple virtual sites are all calculated from the same particles.
*/
void testOverlappingSites() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
NonbondedForce* nonbonded = new NonbondedForce();
system.addForce(nonbonded);
nonbonded->addParticle(1.0, 0.0, 0.0);
nonbonded->addParticle(-0.5, 0.0, 0.0);
nonbonded->addParticle(-0.5, 0.0, 0.0);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(10, 0, 0));
positions.push_back(Vec3(0, 10, 0));
for (int i = 0; i < 20; i++) {
system.addParticle(0.0);
double u = 0.1*((i+1)%4);
double v = 0.05*i;
system.setVirtualSite(3+i, new ThreeParticleAverageSite(0, 1, 2, u, v, 1-u-v));
nonbonded->addParticle(i%2 == 0 ? -1.0 : 1.0, 0.0, 0.0);
positions.push_back(Vec3());
}
VerletIntegrator i1(0.002);
VerletIntegrator i2(0.002);
Context c1(system, i1, Platform::getPlatformByName("Reference"));
Context c2(system, i2, platform);
c1.setPositions(positions);
c2.setPositions(positions);
c1.applyConstraints(0.0001);
c2.applyConstraints(0.0001);
State s1 = c1.getState(State::Positions | State::Forces);
State s2 = c2.getState(State::Positions | State::Forces);
for (int i = 0; i < system.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getPositions()[i], s2.getPositions()[i], 1e-5);
for (int i = 0; i < 3; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
void runPlatformTests();
int main(int argc, char* argv[]) {
try {
initializeTests(argc, argv);
testMasslessParticle();
testTwoParticleAverage();
testThreeParticleAverage();
testOutOfPlane();
testLocalCoordinates();
testConservationLaws();
testOverlappingSites();
runPlatformTests();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
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