Skip to content

GitLab

Commit 2f553a66 authored by Peter Eastman's avatar Peter Eastman
Browse files

Continuing to refactor tests

parent 4ed151e1
Hide whitespace changes
Inline Side-by-side
Showing with 180 additions and 5952 deletions
platforms/cpu/tests/CpuTests.h
View file @ 2f553a66
......@@ -30,8 +30,14 @@
* -------------------------------------------------------------------------- */
#include "CpuPlatform.h"
#include <cstdlib>
#include <iostream>
OpenMM::CpuPlatform platform;
void initializeTests(int argc, char* argv[]) {
if (!OpenMM::CpuPlatform::isProcessorSupported()) {
std::cout << "CPU is not supported. Exiting." << std::endl;
exit(0);
}
}
platforms/cpu/tests/TestCpuCheckpoints.cpp 0 → 100644
View file @ 2f553a66
/* -------------------------------------------------------------------------- *
* 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 "CpuTests.h"
#include "TestCheckpoints.h"
void testCheckpoint() {
const int numParticles = 100;
const double boxSize = 5.0;
const double temperature = 200.0;
System system;
system.addForce(new AndersenThermostat(0.0, 100.0));
NonbondedForce* nonbonded = new NonbondedForce();
system.addForce(nonbonded);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
nonbonded->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.2, 0.1);
bool clash;
do {
clash = false;
positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
for (int j = 0; j < i; j++) {
Vec3 delta = positions[i]-positions[j];
if (sqrt(delta.dot(delta)) < 0.1)
clash = true;
}
} while (clash);
}
VerletIntegrator integrator(0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
context.setParameter(AndersenThermostat::Temperature(), temperature);
// Run for a little while.
integrator.step(100);
// Record the current state and make a checkpoint.
State s1 = context.getState(State::Positions | State::Velocities | State::Parameters);
stringstream stream1(ios_base::out | ios_base::in | ios_base::binary);
context.createCheckpoint(stream1);
// Continue the simulation for a few more steps and record the state again.
integrator.step(10);
State s2 = context.getState(State::Positions | State::Velocities | State::Parameters);
// Restore from the checkpoint and see if everything gets restored correctly.
context.setPeriodicBoxVectors(Vec3(2*boxSize, 0, 0), Vec3(0, 2*boxSize, 0), Vec3(0, 0, 2*boxSize));
context.setParameter(AndersenThermostat::Temperature(), temperature+10);
context.loadCheckpoint(stream1);
State s3 = context.getState(State::Positions | State::Velocities | State::Parameters);
compareStates(s1, s3);
// Now simulate from there and see if the trajectory is identical.
integrator.step(10);
State s4 = context.getState(State::Positions | State::Velocities | State::Parameters);
compareStates(s2, s4);
}
void runPlatformTests() {
testCheckpoint();
}
platforms/cpu/tests/TestCpuEwald.cpp
View file @ 2f553a66
......@@ -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-2013 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 Ewald summation method CPU implementation of NonbondedForce.
*/
#include "CpuTests.h"
#include "TestEwald.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CpuPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/NonbondedForceImpl.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CpuPlatform platform;
const double TOL = 1e-5;
void testEwaldPME(bool includeExceptions) {
// Use amorphous NaCl system for the tests
const int numParticles = 894;
const double cutoff = 1.2;
const double boxSize = 3.00646;
double tol = 1e-5;
ReferencePlatform reference;
System system;
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(tol);
for (int i = 0; i < numParticles/2; i++)
system.addParticle(22.99);
for (int i = 0; i < numParticles/2; i++)
system.addParticle(35.45);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(1.0, 1.0,0.0);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(-1.0, 1.0,0.0);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
vector<Vec3> positions(numParticles);
#include "nacl_amorph.dat"
if (includeExceptions) {
// Add some exclusions.
for (int i = 0; i < numParticles-1; i++) {
Vec3 delta = positions[i]-positions[i+1];
if (sqrt(delta.dot(delta)) < 0.5*cutoff)
nonbonded->addException(i, i+1, i%2 == 0 ? 0.0 : 0.5, 1.0, 0.0);
}
}
// (1) Check whether the Reference and CPU platforms agree when using Ewald Method
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context cpuContext(system, integrator1, platform);
Context referenceContext(system, integrator2, reference);
cpuContext.setPositions(positions);
referenceContext.setPositions(positions);
State cpuState = cpuContext.getState(State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Forces | State::Energy);
tol = 1e-2;
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(referenceState.getForces()[i], cpuState.getForces()[i], tol);
}
tol = 1e-5;
ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), cpuState.getPotentialEnergy(), tol);
// (2) Check whether Ewald method in CPU is self-consistent
double norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = cpuState.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 5e-3;
double step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = cpuState.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
VerletIntegrator integrator3(0.01);
Context cpuContext2(system, integrator3, platform);
cpuContext2.setPositions(positions);
tol = 1e-2;
State cpuState2 = cpuContext2.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (cpuState2.getPotentialEnergy()-cpuState.getPotentialEnergy())/delta, tol)
// (3) Check whether the Reference and CPU platforms agree when using PME
nonbonded->setNonbondedMethod(NonbondedForce::PME);
cpuContext.reinitialize();
referenceContext.reinitialize();
cpuContext.setPositions(positions);
referenceContext.setPositions(positions);
cpuState = cpuContext.getState(State::Forces | State::Energy);
referenceState = referenceContext.getState(State::Forces | State::Energy);
tol = 1e-2;
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(referenceState.getForces()[i], cpuState.getForces()[i], tol);
}
tol = 1e-5;
ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), cpuState.getPotentialEnergy(), tol);
// (4) Check whether PME method in CPU is self-consistent
norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = cpuState.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = cpuState.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
VerletIntegrator integrator4(0.01);
Context cpuContext3(system, integrator4, platform);
cpuContext3.setPositions(positions);
tol = 1e-2;
State cpuState3 = cpuContext3.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (cpuState3.getPotentialEnergy()-cpuState.getPotentialEnergy())/delta, tol)
}
void testEwald2Ions() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(-1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(TOL);
system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(3.048000,2.764000,3.156000);
positions[1] = Vec3(2.809000,2.888000,2.571000);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
ASSERT_EQUAL_VEC(Vec3( 123.711, -64.1877, 302.716), forces[1], 10*TOL);
ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 0.01/*10*TOL*/);
}
void testTriclinic() {
// Create a triclinic box containing eight particles.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(2.5, 0, 0), Vec3(0.5, 3.0, 0), Vec3(0.7, 0.9, 3.5));
for (int i = 0; i < 8; i++)
system.addParticle(1.0);
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
force->setNonbondedMethod(NonbondedForce::PME);
force->setCutoffDistance(1.0);
force->setPMEParameters(3.45891, 32, 40, 48);
for (int i = 0; i < 4; i++)
force->addParticle(-1, 0.440104, 0.4184); // Cl parameters
for (int i = 0; i < 4; i++)
force->addParticle(1, 0.332840, 0.0115897); // Na parameters
vector<Vec3> positions(8);
positions[0] = Vec3(1.744, 2.788, 3.162);
positions[1] = Vec3(1.048, 0.762, 2.340);
positions[2] = Vec3(2.489, 1.570, 2.817);
positions[3] = Vec3(1.027, 1.893, 3.271);
positions[4] = Vec3(0.937, 0.825, 0.009);
positions[5] = Vec3(2.290, 1.887, 3.352);
positions[6] = Vec3(1.266, 1.111, 2.894);
positions[7] = Vec3(0.933, 1.862, 3.490);
// Compute the forces and energy.
VerletIntegrator integ(0.001);
Context context(system, integ, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
// Compare them to values computed by Gromacs.
double expectedEnergy = -963.370;
vector<Vec3> expectedForce(8);
expectedForce[0] = Vec3(4.25253e+01, -1.23503e+02, 1.22139e+02);
expectedForce[1] = Vec3(9.74752e+01, 1.68213e+02, 1.93169e+02);
expectedForce[2] = Vec3(-1.50348e+02, 1.29165e+02, 3.70435e+02);
expectedForce[3] = Vec3(9.18644e+02, -3.52571e+00, -1.34772e+03);
expectedForce[4] = Vec3(-1.61193e+02, 9.01528e+01, -7.12904e+01);
expectedForce[5] = Vec3(2.82630e+02, 2.78029e+01, -3.72864e+02);
expectedForce[6] = Vec3(-1.47454e+02, -2.14448e+02, -3.55789e+02);
expectedForce[7] = Vec3(-8.82195e+02, -7.39132e+01, 1.46202e+03);
for (int i = 0; i < 8; i++) {
ASSERT_EQUAL_VEC(expectedForce[i], state.getForces()[i], 1e-4);
}
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-4);
}
void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
// Create a cloud of random point charges.
const int numParticles = 51;
const double boxWidth = 5.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
}
force->setNonbondedMethod(method);
// For various values of the cutoff and error tolerance, see if the actual error is reasonable.
for (double cutoff = 1.0; cutoff < boxWidth/2; cutoff *= 1.2) {
force->setCutoffDistance(cutoff);
vector<Vec3> refForces;
double norm = 0.0;
for (double tol = 5e-5; tol < 1e-3; tol *= 2.0) {
force->setEwaldErrorTolerance(tol);
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces);
if (refForces.size() == 0) {
refForces = state.getForces();
for (int i = 0; i < numParticles; i++)
norm += refForces[i].dot(refForces[i]);
norm = sqrt(norm);
}
else {
double diff = 0.0;
for (int i = 0; i < numParticles; i++) {
Vec3 delta = refForces[i]-state.getForces()[i];
diff += delta.dot(delta);
}
diff = sqrt(diff)/norm;
ASSERT(diff < 2*tol);
}
if (method == NonbondedForce::PME) {
// See if the PME parameters were calculated correctly.
double expectedAlpha, actualAlpha;
int expectedSize[3], actualSize[3];
NonbondedForceImpl::calcPMEParameters(system, *force, expectedAlpha, expectedSize[0], expectedSize[1], expectedSize[2]);
force->getPMEParametersInContext(context, actualAlpha, actualSize[0], actualSize[1], actualSize[2]);
ASSERT_EQUAL_TOL(expectedAlpha, actualAlpha, 1e-5);
for (int i = 0; i < 3; i++) {
ASSERT(actualSize[i] >= expectedSize[i]);
ASSERT(actualSize[i] < expectedSize[i]+10);
}
}
}
}
}
int main(int argc, char* argv[]) {
try {
if (!CpuPlatform::isProcessorSupported()) {
cout << "CPU is not supported. Exiting." << endl;
return 0;
}
testEwaldPME(false);
testEwaldPME(true);
// testEwald2Ions();
testTriclinic();
testErrorTolerance(NonbondedForce::Ewald);
testErrorTolerance(NonbondedForce::PME);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cpu/tests/TestCpuGBSAOBCForce.cpp
View file @ 2f553a66
......@@ -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-2014 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,244 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CPU implementation of GBSAOBCForce.
*/
#include "CpuTests.h"
#include "TestGBSAOBCForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CpuPlatform.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/NonbondedForce.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testSingleParticle() {
CpuPlatform platform;
System system;
system.addParticle(2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* forceField = new GBSAOBCForce();
forceField->addParticle(0.5, 0.15, 1);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double bornRadius = 0.15-0.009; // dielectric offset
double eps0 = EPSILON0;
double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric())/bornRadius;
double extendedRadius = 0.15+0.14; // probe radius
double nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
// Change the parameters and see if it is still correct.
forceField->setParticleParameters(0, 0.4, 0.25, 1);
forceField->updateParametersInContext(context);
state = context.getState(State::Energy);
bornRadius = 0.25-0.009; // dielectric offset
bornEnergy = (-0.4*0.4/(8*PI_M*eps0))*(1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric())/bornRadius;
extendedRadius = 0.25+0.14;
nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.25/bornRadius, 6.0);
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
}
void testGlobalSettings() {
CpuPlatform platform;
System system;
system.addParticle(2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* forceField = new GBSAOBCForce();
forceField->addParticle(0.5, 0.15, 1);
const double soluteDielectric = 2.1;
const double solventDielectric = 35.0;
const double surfaceAreaEnergy = 0.75;
forceField->setSoluteDielectric(soluteDielectric);
forceField->setSolventDielectric(solventDielectric);
forceField->setSurfaceAreaEnergy(surfaceAreaEnergy);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double bornRadius = 0.15-0.009; // dielectric offset
double eps0 = EPSILON0;
double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/soluteDielectric-1.0/solventDielectric)/bornRadius;
double extendedRadius = 0.15+0.14; // probe radius
double nonpolarEnergy = 4*PI_M*surfaceAreaEnergy*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
}
void testCutoffAndPeriodic() {
CpuPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* gbsa = new GBSAOBCForce();
NonbondedForce* nonbonded = new NonbondedForce();
gbsa->addParticle(-1, 0.15, 1);
nonbonded->addParticle(-1, 1, 0);
gbsa->addParticle(1, 0.15, 1);
nonbonded->addParticle(1, 1, 0);
const double cutoffDistance = 3.0;
const double boxSize = 10.0;
nonbonded->setCutoffDistance(cutoffDistance);
gbsa->setCutoffDistance(cutoffDistance);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(gbsa);
system.addForce(nonbonded);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
// Calculate the forces for both cutoff and periodic with two different atom positions.
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
Context context(system, integrator, platform);
context.setPositions(positions);
State state1 = context.getState(State::Forces);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
context.reinitialize();
context.setPositions(positions);
State state2 = context.getState(State::Forces);
positions[1][0]+= boxSize;
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
context.reinitialize();
context.setPositions(positions);
State state3 = context.getState(State::Forces);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
context.reinitialize();
context.setPositions(positions);
State state4 = context.getState(State::Forces);
// All forces should be identical, exception state3 which should be zero.
ASSERT_EQUAL_VEC(state1.getForces()[0], state2.getForces()[0], 0.01);
ASSERT_EQUAL_VEC(state1.getForces()[1], state2.getForces()[1], 0.01);
ASSERT_EQUAL_VEC(state1.getForces()[0], state4.getForces()[0], 0.01);
ASSERT_EQUAL_VEC(state1.getForces()[1], state4.getForces()[1], 0.01);
ASSERT_EQUAL_VEC(state3.getForces()[0], Vec3(0, 0, 0), 0.01);
ASSERT_EQUAL_VEC(state3.getForces()[1], Vec3(0, 0, 0), 0.01);
}
void testForce(int numParticles, NonbondedForce::NonbondedMethod method, GBSAOBCForce::NonbondedMethod method2) {
CpuPlatform platform;
ReferencePlatform reference;
System system;
GBSAOBCForce* gbsa = new GBSAOBCForce();
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
double charge = i%2 == 0 ? -1 : 1;
gbsa->addParticle(charge, 0.15, 1);
nonbonded->addParticle(charge, 1, 0);
}
nonbonded->setNonbondedMethod(method);
gbsa->setNonbondedMethod(method2);
nonbonded->setCutoffDistance(3.0);
gbsa->setCutoffDistance(3.0);
int grid = (int) floor(0.5+pow(numParticles, 1.0/3.0));
if (method == NonbondedForce::CutoffPeriodic) {
double boxSize = (grid+1)*1.1;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
}
system.addForce(gbsa);
system.addForce(nonbonded);
LangevinIntegrator integrator1(0, 0.1, 0.01);
LangevinIntegrator integrator2(0, 0.1, 0.01);
Context context(system, integrator1, platform);
Context refContext(system, integrator2, reference);
// Set random (but uniformly distributed) positions for all the particles.
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < grid; i++)
for (int j = 0; j < grid; j++)
for (int k = 0; k < grid; k++)
positions[i*grid*grid+j*grid+k] = Vec3(i*1.1, j*1.1, k*1.1);
for (int i = 0; i < numParticles; ++i)
positions[i] = positions[i] + Vec3(0.5*genrand_real2(sfmt), 0.5*genrand_real2(sfmt), 0.5*genrand_real2(sfmt));
context.setPositions(positions);
refContext.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
State refState = refContext.getState(State::Forces | State::Energy);
// Make sure the CPU and Reference platforms agree.
double norm = 0.0;
double diff = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = state.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
Vec3 delta = f-refState.getForces()[i];
diff += delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
}
norm = std::sqrt(norm);
diff = std::sqrt(diff);
ASSERT_EQUAL_TOL(0.0, diff, 0.001*norm);
ASSERT_EQUAL_TOL(state.getPotentialEnergy(), refState.getPotentialEnergy(), 1e-3);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
// (This doesn't work with cutoffs, since the energy changes discontinuously at the cutoff distance.)
if (method == NonbondedForce::NoCutoff)
{
const double delta = 0.3;
double step = 0.5*delta/norm;
vector<Vec3> positions2(numParticles), positions3(numParticles);
for (int i = 0; i < numParticles; ++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);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/delta, 1e-2)
}
}
int main() {
try {
if (!CpuPlatform::isProcessorSupported()) {
cout << "CPU is not supported. Exiting." << endl;
return 0;
}
testSingleParticle();
testGlobalSettings();
testCutoffAndPeriodic();
for (int i = 5; i < 11; i++) {
testForce(i*i*i, NonbondedForce::NoCutoff, GBSAOBCForce::NoCutoff);
testForce(i*i*i, NonbondedForce::CutoffNonPeriodic, GBSAOBCForce::CutoffNonPeriodic);
testForce(i*i*i, NonbondedForce::CutoffPeriodic, GBSAOBCForce::CutoffPeriodic);
}
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cpu/tests/TestCpuLangevinIntegrator.cpp
View file @ 2f553a66
......@@ -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-2013 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,255 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of LangevinIntegrator.
*/
#include "CpuTests.h"
#include "TestLangevinIntegrator.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CpuPlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.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() {
CpuPlatform platform;
System system;
system.addParticle(2.0);
system.addParticle(2.0);
LangevinIntegrator integrator(0, 0.1, 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 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);
}
// Not set the friction to a tiny value and see if it conserves energy.
integrator.setFriction(5e-5);
context.setPositions(positions);
State state = context.getState(State::Energy);
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.01);
integrator.step(1);
}
}
void testTemperature() {
const int numParticles = 8;
const double temp = 100.0;
CpuPlatform platform;
System system;
LangevinIntegrator integrator(temp, 2.0, 0.01);
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);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 10000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
}
ke /= 10000;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 6/std::sqrt(10000.0));
}
void testConstraints() {
const int numParticles = 8;
const double temp = 100.0;
CpuPlatform platform;
System system;
LangevinIntegrator integrator(temp, 2.0, 0.01);
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);
}
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() {
CpuPlatform platform;
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);
LangevinIntegrator 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;
const double collisionFreq = 10.0;
CpuPlatform platform;
System system;
LangevinIntegrator integrator(temp, 2.0, 0.01);
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_EQUAL_TOL(state1.getPositions()[i][j], state2.getPositions()[i][j], 1e-5);
ASSERT_EQUAL_TOL(state3.getPositions()[i][j], state4.getPositions()[i][j], 1e-5);
ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
}
}
}
int main() {
try {
if (!CpuPlatform::isProcessorSupported()) {
cout << "CPU is not supported. Exiting." << endl;
return 0;
}
testSingleBond();
testTemperature();
testConstraints();
testConstrainedMasslessParticles();
testRandomSeed();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cpu/tests/TestCpuNonbondedForce.cpp
View file @ 2f553a66
......@@ -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-2015 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,684 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests all the different force terms in the CUDA implementation of NonbondedForce.
*/
#include "CpuTests.h"
#include "TestNonbondedForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CpuPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/internal/ContextImpl.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CpuPlatform platform;
const double TOL = 1e-5;
void testCoulomb() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->addParticle(0.5, 1, 0);
forceField->addParticle(-1.5, 1, 0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double force = ONE_4PI_EPS0*(-0.75)/4.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(-0.75)/2.0, state.getPotentialEnergy(), TOL);
}
void testLJ() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->addParticle(0, 1.2, 1);
forceField->addParticle(0, 1.4, 2);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double x = 1.3/2.0;
double eps = SQRT_TWO;
double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/2.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0)), state.getPotentialEnergy(), TOL);
}
void testExclusionsAnd14() {
System system;
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < 5; ++i) {
system.addParticle(1.0);
nonbonded->addParticle(0, 1.5, 0);
}
vector<pair<int, int> > bonds;
bonds.push_back(pair<int, int>(0, 1));
bonds.push_back(pair<int, int>(1, 2));
bonds.push_back(pair<int, int>(2, 3));
bonds.push_back(pair<int, int>(3, 4));
nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
int first14, second14;
for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon;
nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
first14 = i;
if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
second14 = i;
}
system.addForce(nonbonded);
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
for (int i = 1; i < 5; ++i) {
// Test LJ forces
vector<Vec3> positions(5);
const double r = 1.0;
for (int j = 0; j < 5; ++j) {
nonbonded->setParticleParameters(j, 0, 1.5, 0);
positions[j] = Vec3(0, j, 0);
}
nonbonded->setParticleParameters(0, 0, 1.5, 1);
nonbonded->setParticleParameters(i, 0, 1.5, 1);
nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
positions[i] = Vec3(r, 0, 0);
context.reinitialize();
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double x = 1.5/r;
double eps = 1.0;
double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
double energy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
if (i == 3) {
force *= 0.5;
energy *= 0.5;
}
if (i < 3) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
// Test Coulomb forces
nonbonded->setParticleParameters(0, 2, 1.5, 0);
nonbonded->setParticleParameters(i, 2, 1.5, 0);
nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces();
force = ONE_4PI_EPS0*4/(r*r);
energy = ONE_4PI_EPS0*4/r;
if (i == 3) {
force /= 1.2;
energy /= 1.2;
}
if (i < 3) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
}
void testCutoff() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->addParticle(1.0, 1, 0);
forceField->addParticle(1.0, 1, 0);
forceField->addParticle(1.0, 1, 0);
forceField->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
const double cutoff = 2.9;
forceField->setCutoffDistance(cutoff);
const double eps = 50.0;
forceField->setReactionFieldDielectric(eps);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(0, 2, 0);
positions[2] = Vec3(0, 3, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force1 = ONE_4PI_EPS0*(1.0)*(0.25-2.0*krf*2.0);
const double force2 = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL);
const double energy1 = ONE_4PI_EPS0*(1.0)*(0.5+krf*4.0-crf);
const double energy2 = ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf);
ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL);
}
void testCutoff14() {
System system;
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
for (int i = 0; i < 5; ++i) {
system.addParticle(1.0);
nonbonded->addParticle(0, 1.5, 0);
}
const double cutoff = 3.5;
nonbonded->setCutoffDistance(cutoff);
const double eps = 30.0;
nonbonded->setReactionFieldDielectric(eps);
vector<pair<int, int> > bonds;
bonds.push_back(pair<int, int>(0, 1));
bonds.push_back(pair<int, int>(1, 2));
bonds.push_back(pair<int, int>(2, 3));
bonds.push_back(pair<int, int>(3, 4));
nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
int first14, second14;
for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon;
nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
first14 = i;
if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
second14 = i;
}
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(5);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(2, 0, 0);
positions[3] = Vec3(3, 0, 0);
positions[4] = Vec3(4, 0, 0);
for (int i = 1; i < 5; ++i) {
// Test LJ forces
nonbonded->setParticleParameters(0, 0, 1.5, 1);
for (int j = 1; j < 5; ++j)
nonbonded->setParticleParameters(j, 0, 1.5, 0);
nonbonded->setParticleParameters(i, 0, 1.5, 1);
nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
context.reinitialize();
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double r = positions[i][0];
double x = 1.5/r;
double e = 1.0;
double force = 4.0*e*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
double energy = 4.0*e*(std::pow(x, 12.0)-std::pow(x, 6.0));
if (i == 3) {
force *= 0.5;
energy *= 0.5;
}
if (i < 3 || r > cutoff) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
// Test Coulomb forces
const double q = 0.7;
nonbonded->setParticleParameters(0, q, 1.5, 0);
nonbonded->setParticleParameters(i, q, 1.5, 0);
nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces();
force = ONE_4PI_EPS0*q*q/(r*r);
energy = ONE_4PI_EPS0*q*q/r;
if (i == 3) {
force /= 1.2;
energy /= 1.2;
}
if (i < 3 || r > cutoff) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
}
void testPeriodic() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addException(0, 1, 0.0, 1.0, 0.0);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
positions[2] = Vec3(3, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
const double eps = 78.3;
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
}
void testTriclinic() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
Vec3 a(3.1, 0, 0);
Vec3 b(0.4, 3.5, 0);
Vec3 c(-0.1, -0.5, 4.0);
system.setDefaultPeriodicBoxVectors(a, b, c);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 1.5;
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const double eps = 78.3;
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
for (int iteration = 0; iteration < 50; iteration++) {
// Generate random positions for the two particles.
positions[0] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
positions[1] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
context.setPositions(positions);
// Loop over all possible periodic copies and find the nearest one.
Vec3 delta;
double distance2 = 100.0;
for (int i = -1; i < 2; i++)
for (int j = -1; j < 2; j++)
for (int k = -1; k < 2; k++) {
Vec3 d = positions[1]-positions[0]+a*i+b*j+c*k;
if (d.dot(d) < distance2) {
delta = d;
distance2 = d.dot(d);
}
}
double distance = sqrt(distance2);
// See if the force and energy are correct.
State state = context.getState(State::Forces | State::Energy);
if (distance >= cutoff) {
ASSERT_EQUAL(0.0, state.getPotentialEnergy());
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getForces()[0], 0);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getForces()[1], 0);
}
else {
const Vec3 force = delta*ONE_4PI_EPS0*(-1.0/(distance*distance*distance)+2.0*krf);
ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(1.0/distance+krf*distance*distance-crf), state.getPotentialEnergy(), 1e-4);
ASSERT_EQUAL_VEC(force, state.getForces()[0], 2e-5);
ASSERT_EQUAL_VEC(-force, state.getForces()[1], 2e-5);
}
}
}
void testLargeSystem() {
const int numMolecules = 600;
const int numParticles = numMolecules*2;
const double cutoff = 2.0;
const double boxSize = 20.0;
const double tol = 2e-3;
ReferencePlatform reference;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
NonbondedForce* nonbonded = new NonbondedForce();
HarmonicBondForce* bonds = new HarmonicBondForce();
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numMolecules; i++) {
if (i < numMolecules/2) {
nonbonded->addParticle(-1.0, 0.2, 0.1);
nonbonded->addParticle(1.0, 0.1, 0.1);
}
else {
nonbonded->addParticle(-1.0, 0.2, 0.2);
nonbonded->addParticle(1.0, 0.1, 0.2);
}
positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
velocities[2*i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
velocities[2*i+1] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
bonds->addBond(2*i, 2*i+1, 1.0, 0.1);
nonbonded->addException(2*i, 2*i+1, 0.0, 0.15, 0.0);
}
// Try with cutoffs but not periodic boundary conditions, and make sure the cl and Reference
// platforms agree.
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
system.addForce(bonds);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context cpuContext(system, integrator1, platform);
Context referenceContext(system, integrator2, reference);
cpuContext.setPositions(positions);
cpuContext.setVelocities(velocities);
referenceContext.setPositions(positions);
referenceContext.setVelocities(velocities);
State cpuState = cpuContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(cpuState.getPositions()[i], referenceState.getPositions()[i], tol);
ASSERT_EQUAL_VEC(cpuState.getVelocities()[i], referenceState.getVelocities()[i], tol);
ASSERT_EQUAL_VEC(cpuState.getForces()[i], referenceState.getForces()[i], tol);
}
ASSERT_EQUAL_TOL(cpuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
// Now do the same thing with periodic boundary conditions.
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
cpuContext.reinitialize();
referenceContext.reinitialize();
cpuContext.setPositions(positions);
cpuContext.setVelocities(velocities);
referenceContext.setPositions(positions);
referenceContext.setVelocities(velocities);
cpuState = cpuContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++) {
double dx = cpuState.getPositions()[i][0]-referenceState.getPositions()[i][0];
double dy = cpuState.getPositions()[i][1]-referenceState.getPositions()[i][1];
double dz = cpuState.getPositions()[i][2]-referenceState.getPositions()[i][2];
ASSERT_EQUAL_TOL(fmod(cpuState.getPositions()[i][0]-referenceState.getPositions()[i][0], boxSize), 0, tol);
ASSERT_EQUAL_TOL(fmod(cpuState.getPositions()[i][1]-referenceState.getPositions()[i][1], boxSize), 0, tol);
ASSERT_EQUAL_TOL(fmod(cpuState.getPositions()[i][2]-referenceState.getPositions()[i][2], boxSize), 0, tol);
ASSERT_EQUAL_VEC(cpuState.getVelocities()[i], referenceState.getVelocities()[i], tol);
ASSERT_EQUAL_VEC(cpuState.getForces()[i], referenceState.getForces()[i], tol);
}
ASSERT_EQUAL_TOL(cpuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
}
void testDispersionCorrection() {
// Create a box full of identical particles.
int gridSize = 5;
int numParticles = gridSize*gridSize*gridSize;
double boxSize = gridSize*0.7;
double cutoff = boxSize/3;
System system;
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions(numParticles);
int index = 0;
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
system.addParticle(1.0);
nonbonded->addParticle(0, 1.1, 0.5);
positions[index] = Vec3(i*boxSize/gridSize, j*boxSize/gridSize, k*boxSize/gridSize);
index++;
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
// See if the correction has the correct value.
Context context(system, integrator, platform);
context.setPositions(positions);
double energy1 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseDispersionCorrection(false);
context.reinitialize();
context.setPositions(positions);
double energy2 = context.getState(State::Energy).getPotentialEnergy();
double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
// Now modify half the particles to be different, and see if it is still correct.
int numType2 = 0;
for (int i = 0; i < numParticles; i += 2) {
nonbonded->setParticleParameters(i, 0, 1, 1);
numType2++;
}
int numType1 = numParticles-numType2;
nonbonded->updateParametersInContext(context);
energy2 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseDispersionCorrection(true);
context.reinitialize();
context.setPositions(positions);
energy1 = context.getState(State::Energy).getPotentialEnergy();
term1 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
term2 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
term1 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 12)/pow(cutoff, 9))/9;
term2 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 6)/pow(cutoff, 3))/3;
double combinedSigma = 0.5*(1+1.1);
double combinedEpsilon = sqrt(1*0.5);
term1 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 12)/pow(cutoff, 9))/9;
term2 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 6)/pow(cutoff, 3))/3;
term1 /= (numParticles*(numParticles+1))/2;
term2 /= (numParticles*(numParticles+1))/2;
expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
}
void testChangingParameters() {
const int numMolecules = 600;
const int numParticles = numMolecules*2;
const double cutoff = 2.0;
const double boxSize = 20.0;
const double tol = 2e-3;
ReferencePlatform reference;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numMolecules; i++) {
if (i < numMolecules/2) {
nonbonded->addParticle(-1.0, 0.2, 0.1);
nonbonded->addParticle(1.0, 0.1, 0.1);
}
else {
nonbonded->addParticle(-1.0, 0.2, 0.2);
nonbonded->addParticle(1.0, 0.1, 0.2);
}
positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
system.addConstraint(2*i, 2*i+1, 1.0);
nonbonded->addException(2*i, 2*i+1, 0.0, 0.15, 0.0);
}
nonbonded->setNonbondedMethod(NonbondedForce::PME);
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
// See if Reference and CPU give the same forces and energies.
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context cpuContext(system, integrator1, platform);
Context referenceContext(system, integrator2, reference);
cpuContext.setPositions(positions);
referenceContext.setPositions(positions);
State cpuState = cpuContext.getState(State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(cpuState.getForces()[i], referenceState.getForces()[i], tol);
ASSERT_EQUAL_TOL(cpuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
// Now modify parameters and see if they still agree.
for (int i = 0; i < numParticles; i += 5) {
double charge, sigma, epsilon;
nonbonded->getParticleParameters(i, charge, sigma, epsilon);
nonbonded->setParticleParameters(i, 1.5*charge, 1.1*sigma, 1.7*epsilon);
}
nonbonded->updateParametersInContext(cpuContext);
nonbonded->updateParametersInContext(referenceContext);
cpuState = cpuContext.getState(State::Forces | State::Energy);
referenceState = referenceContext.getState(State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(cpuState.getForces()[i], referenceState.getForces()[i], tol);
ASSERT_EQUAL_TOL(cpuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
}
void testSwitchingFunction(NonbondedForce::NonbondedMethod method) {
System system;
system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(0, 1.2, 1);
nonbonded->addParticle(0, 1.4, 2);
nonbonded->setNonbondedMethod(method);
nonbonded->setCutoffDistance(2.0);
nonbonded->setUseSwitchingFunction(true);
nonbonded->setSwitchingDistance(1.5);
nonbonded->setUseDispersionCorrection(false);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
double eps = SQRT_TWO;
// Compute the interaction at various distances.
for (double r = 1.0; r < 2.5; r += 0.1) {
positions[1] = Vec3(r, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
// See if the energy is correct.
double x = 1.3/r;
double expectedEnergy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
double switchValue;
if (r <= 1.5)
switchValue = 1;
else if (r >= 2.0)
switchValue = 0;
else {
double t = (r-1.5)/0.5;
switchValue = 1+t*t*t*(-10+t*(15-t*6));
}
ASSERT_EQUAL_TOL(switchValue*expectedEnergy, state.getPotentialEnergy(), TOL);
// See if the force is the gradient of the energy.
double delta = 1e-3;
positions[1] = Vec3(r-delta, 0, 0);
context.setPositions(positions);
double e1 = context.getState(State::Energy).getPotentialEnergy();
positions[1] = Vec3(r+delta, 0, 0);
context.setPositions(positions);
double e2 = context.getState(State::Energy).getPotentialEnergy();
ASSERT_EQUAL_TOL((e2-e1)/(2*delta), state.getForces()[0][0], 1e-3);
}
}
int main(int argc, char* argv[]) {
try {
if (!CpuPlatform::isProcessorSupported()) {
cout << "CPU is not supported. Exiting." << endl;
return 0;
}
testCoulomb();
testLJ();
testExclusionsAnd14();
testCutoff();
testCutoff14();
testPeriodic();
testTriclinic();
testLargeSystem();
testDispersionCorrection();
testChangingParameters();
testSwitchingFunction(NonbondedForce::CutoffNonPeriodic);
testSwitchingFunction(NonbondedForce::PME);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cpu/tests/TestCpuSettle.cpp
View file @ 2f553a66
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
......@@ -7,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-2013 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -30,91 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CPU implementation of the SETTLE algorithm.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CpuPlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
void testConstraints() {
const int numMolecules = 10;
const int numParticles = numMolecules*3;
const int numConstraints = numMolecules*3;
const double temp = 100.0;
CpuPlatform platform;
System system;
LangevinIntegrator integrator(temp, 2.0, 0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numMolecules; ++i) {
system.addParticle(16.0);
system.addParticle(1.0);
system.addParticle(1.0);
forceField->addParticle(-0.82, 0.317, 0.65);
forceField->addParticle(0.41, 1.0, 0.0);
forceField->addParticle(0.41, 1.0, 0.0);
system.addConstraint(i*3, i*3+1, 0.1);
system.addConstraint(i*3, i*3+2, 0.1);
system.addConstraint(i*3+1, i*3+2, 0.163);
}
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 < numMolecules; ++i) {
positions[i*3] = Vec3((i%4)*0.4, (i/4)*0.4, 0);
positions[i*3+1] = positions[i*3]+Vec3(0.1, 0, 0);
positions[i*3+2] = positions[i*3]+Vec3(-0.03333, 0.09428, 0);
velocities[i*3] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
velocities[i*3+1] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
velocities[i*3+2] = 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) {
integrator.step(1);
State state = context.getState(State::Positions | 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-5);
}
}
}
#include "CpuTests.h"
#include "TestSettle.h"
int main(int argc, char* argv[]) {
try {
if (!CpuPlatform::isProcessorSupported()) {
cout << "CPU is not supported. Exiting." << endl;
return 0;
}
testConstraints();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cuda/tests/TestCudaCheckpoints.cpp
View file @ 2f553a66
......@@ -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-2013 Stanford University and the Authors. *
* Portions copyright (c) 2012-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,45 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests creating and loading checkpoints with the CUDA platform.
*/
#include "CudaPlatform.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/AndersenThermostat.h"
#include "openmm/Context.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <sstream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
CudaPlatform platform;
void compareStates(State& s1, State& s2) {
ASSERT_EQUAL_TOL(s1.getTime(), s2.getTime(), TOL);
int numParticles = s1.getPositions().size();
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(s1.getPositions()[i], s2.getPositions()[i], TOL);
ASSERT_EQUAL_VEC(s1.getVelocities()[i], s2.getVelocities()[i], TOL);
Vec3 a1, b1, c1, a2, b2, c2;
s1.getPeriodicBoxVectors(a1, b1, c1);
s2.getPeriodicBoxVectors(a2, b2, c2);
ASSERT_EQUAL_VEC(a1, a2, TOL);
ASSERT_EQUAL_VEC(b1, b2, TOL);
ASSERT_EQUAL_VEC(c1, c2, TOL);
for (map<string, double>::const_iterator iter = s1.getParameters().begin(); iter != s1.getParameters().end(); ++iter)
ASSERT_EQUAL(iter->second, (*s2.getParameters().find(iter->first)).second);
}
}
#include "CudaTests.h"
#include "TestCheckpoints.h"
void testCheckpoint() {
const int numParticles = 100;
......@@ -159,71 +122,6 @@ void testCheckpoint() {
compareStates(s6, s8);
}
void testSetState() {
const int numParticles = 10;
const double boxSize = 3.0;
const double temperature = 200.0;
System system;
system.addForce(new AndersenThermostat(0.0, 100.0));
NonbondedForce* nonbonded = new NonbondedForce();
system.addForce(nonbonded);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
nonbonded->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.2, 0.1);
positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
}
VerletIntegrator integrator(0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
context.setParameter(AndersenThermostat::Temperature(), temperature);
// Run for a little while.
integrator.step(100);
// Record the current state.
State s1 = context.getState(State::Positions | State::Velocities | State::Parameters);
// Continue the simulation for a few more steps and record a partial state.
integrator.step(10);
State s2 = context.getState(State::Positions);
// Restore the original state and see if everything gets restored correctly.
context.setPeriodicBoxVectors(Vec3(2*boxSize, 0, 0), Vec3(0, 2*boxSize, 0), Vec3(0, 0, 2*boxSize));
context.setParameter(AndersenThermostat::Temperature(), temperature+10);
context.setState(s1);
State s3 = context.getState(State::Positions | State::Velocities | State::Parameters);
compareStates(s1, s3);
// Set the partial state and see if the correct things were set.
context.setState(s2);
State s4 = context.getState(State::Positions | State::Velocities | State::Parameters);
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(s2.getPositions()[i], s4.getPositions()[i], TOL);
ASSERT_EQUAL_VEC(s1.getVelocities()[i], s4.getVelocities()[i], TOL);
}
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testCheckpoint();
testSetState();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
testCheckpoint();
}
platforms/cuda/tests/TestCudaEwald.cpp
View file @ 2f553a66
......@@ -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-2014 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,358 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the Ewald summation method CUDA implementation of NonbondedForce.
*/
#include "CudaTests.h"
#include "TestEwald.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/NonbondedForceImpl.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 testEwaldPME(bool includeExceptions) {
// Use amorphous NaCl system for the tests
const int numParticles = 894;
const double cutoff = 1.2;
const double boxSize = 3.00646;
double tol = 1e-5;
ReferencePlatform reference;
System system;
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(tol);
for (int i = 0; i < numParticles/2; i++)
system.addParticle(22.99);
for (int i = 0; i < numParticles/2; i++)
system.addParticle(35.45);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(1.0, 1.0,0.0);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(-1.0, 1.0,0.0);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
vector<Vec3> positions(numParticles);
#include "nacl_amorph.dat"
if (includeExceptions) {
// Add some exclusions.
for (int i = 0; i < numParticles-1; i++) {
Vec3 delta = positions[i]-positions[i+1];
if (sqrt(delta.dot(delta)) < 0.5*cutoff)
nonbonded->addException(i, i+1, i%2 == 0 ? 0.0 : 0.5, 1.0, 0.0);
}
}
// (1) Check whether the Reference and CUDA platforms agree when using Ewald Method
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context cuContext(system, integrator1, platform);
Context referenceContext(system, integrator2, reference);
cuContext.setPositions(positions);
referenceContext.setPositions(positions);
State cuState = cuContext.getState(State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Forces | State::Energy);
tol = 1e-2;
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(referenceState.getForces()[i], cuState.getForces()[i], tol);
}
tol = 1e-5;
ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), cuState.getPotentialEnergy(), tol);
// (2) Check whether Ewald method in CUDA is self-consistent
double norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = cuState.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 5e-3;
double step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = cuState.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
VerletIntegrator integrator3(0.01);
Context cuContext2(system, integrator3, platform);
cuContext2.setPositions(positions);
tol = 1e-2;
State cuState2 = cuContext2.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (cuState2.getPotentialEnergy()-cuState.getPotentialEnergy())/delta, tol)
// (3) Check whether the Reference and CUDA platforms agree when using PME
nonbonded->setNonbondedMethod(NonbondedForce::PME);
cuContext.reinitialize();
referenceContext.reinitialize();
cuContext.setPositions(positions);
referenceContext.setPositions(positions);
cuState = cuContext.getState(State::Forces | State::Energy);
referenceState = referenceContext.getState(State::Forces | State::Energy);
tol = 1e-2;
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(referenceState.getForces()[i], cuState.getForces()[i], tol);
}
tol = 1e-5;
ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), cuState.getPotentialEnergy(), tol);
// (4) Check whether PME method in CUDA is self-consistent
norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = cuState.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = cuState.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
VerletIntegrator integrator4(0.01);
Context cuContext3(system, integrator4, platform);
cuContext3.setPositions(positions);
tol = 1e-2;
State cuState3 = cuContext3.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (cuState3.getPotentialEnergy()-cuState.getPotentialEnergy())/delta, tol)
}
void testEwald2Ions() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(-1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(TOL);
system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(3.048000,2.764000,3.156000);
positions[1] = Vec3(2.809000,2.888000,2.571000);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
ASSERT_EQUAL_VEC(Vec3( 123.711, -64.1877, 302.716), forces[1], 10*TOL);
ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 0.01/*10*TOL*/);
}
void testTriclinic() {
// Create a triclinic box containing eight particles.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(2.5, 0, 0), Vec3(0.5, 3.0, 0), Vec3(0.7, 0.9, 3.5));
for (int i = 0; i < 8; i++)
system.addParticle(1.0);
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
force->setNonbondedMethod(NonbondedForce::PME);
force->setCutoffDistance(1.0);
force->setPMEParameters(3.45891, 32, 40, 48);
for (int i = 0; i < 4; i++)
force->addParticle(-1, 0.440104, 0.4184); // Cl parameters
for (int i = 0; i < 4; i++)
force->addParticle(1, 0.332840, 0.0115897); // Na parameters
vector<Vec3> positions(8);
positions[0] = Vec3(1.744, 2.788, 3.162);
positions[1] = Vec3(1.048, 0.762, 2.340);
positions[2] = Vec3(2.489, 1.570, 2.817);
positions[3] = Vec3(1.027, 1.893, 3.271);
positions[4] = Vec3(0.937, 0.825, 0.009);
positions[5] = Vec3(2.290, 1.887, 3.352);
positions[6] = Vec3(1.266, 1.111, 2.894);
positions[7] = Vec3(0.933, 1.862, 3.490);
// Compute the forces and energy.
VerletIntegrator integ(0.001);
Context context(system, integ, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
// Compare them to values computed by Gromacs.
double expectedEnergy = -963.370;
vector<Vec3> expectedForce(8);
expectedForce[0] = Vec3(4.25253e+01, -1.23503e+02, 1.22139e+02);
expectedForce[1] = Vec3(9.74752e+01, 1.68213e+02, 1.93169e+02);
expectedForce[2] = Vec3(-1.50348e+02, 1.29165e+02, 3.70435e+02);
expectedForce[3] = Vec3(9.18644e+02, -3.52571e+00, -1.34772e+03);
expectedForce[4] = Vec3(-1.61193e+02, 9.01528e+01, -7.12904e+01);
expectedForce[5] = Vec3(2.82630e+02, 2.78029e+01, -3.72864e+02);
expectedForce[6] = Vec3(-1.47454e+02, -2.14448e+02, -3.55789e+02);
expectedForce[7] = Vec3(-8.82195e+02, -7.39132e+01, 1.46202e+03);
for (int i = 0; i < 8; i++) {
ASSERT_EQUAL_VEC(expectedForce[i], state.getForces()[i], 1e-4);
}
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-4);
}
void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
// Create a cloud of random point charges.
const int numParticles = 51;
const double boxWidth = 5.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
}
force->setNonbondedMethod(method);
// For various values of the cutoff and error tolerance, see if the actual error is reasonable.
for (double cutoff = 1.0; cutoff < boxWidth/2; cutoff *= 1.2) {
force->setCutoffDistance(cutoff);
vector<Vec3> refForces;
double norm = 0.0;
for (double tol = 5e-5; tol < 1e-3; tol *= 2.0) {
force->setEwaldErrorTolerance(tol);
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces);
if (refForces.size() == 0) {
refForces = state.getForces();
for (int i = 0; i < numParticles; i++)
norm += refForces[i].dot(refForces[i]);
norm = sqrt(norm);
}
else {
double diff = 0.0;
for (int i = 0; i < numParticles; i++) {
Vec3 delta = refForces[i]-state.getForces()[i];
diff += delta.dot(delta);
}
diff = sqrt(diff)/norm;
ASSERT(diff < 2*tol);
}
if (method == NonbondedForce::PME) {
// See if the PME parameters were calculated correctly.
double expectedAlpha, actualAlpha;
int expectedSize[3], actualSize[3];
NonbondedForceImpl::calcPMEParameters(system, *force, expectedAlpha, expectedSize[0], expectedSize[1], expectedSize[2]);
force->getPMEParametersInContext(context, actualAlpha, actualSize[0], actualSize[1], actualSize[2]);
ASSERT_EQUAL_TOL(expectedAlpha, actualAlpha, 1e-5);
for (int i = 0; i < 3; i++) {
ASSERT(actualSize[i] >= expectedSize[i]);
ASSERT(actualSize[i] < expectedSize[i]+10);
}
}
}
}
}
void testPMEParameters() {
// Create a cloud of random point charges.
const int numParticles = 51;
const double boxWidth = 4.7;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
}
force->setNonbondedMethod(NonbondedForce::PME);
// Compute the energy with an error tolerance of 1e-3.
force->setEwaldErrorTolerance(1e-3);
VerletIntegrator integrator1(0.01);
Context context1(system, integrator1, platform);
context1.setPositions(positions);
double energy1 = context1.getState(State::Energy).getPotentialEnergy();
// Try again with an error tolerance of 1e-4.
force->setEwaldErrorTolerance(1e-4);
VerletIntegrator integrator2(0.01);
Context context2(system, integrator2, platform);
context2.setPositions(positions);
double energy2 = context2.getState(State::Energy).getPotentialEnergy();
// Now explicitly set the parameters. These should match the values that were
// used for tolerance 1e-3.
force->setPMEParameters(2.49291157051793, 32, 32, 32);
VerletIntegrator integrator3(0.01);
Context context3(system, integrator3, platform);
context3.setPositions(positions);
double energy3 = context3.getState(State::Energy).getPotentialEnergy();
ASSERT_EQUAL_TOL(energy1, energy3, 1e-6);
ASSERT(fabs((energy1-energy2)/energy1) > 1e-5);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testEwaldPME(false);
testEwaldPME(true);
// testEwald2Ions();
testTriclinic();
testErrorTolerance(NonbondedForce::Ewald);
testErrorTolerance(NonbondedForce::PME);
testPMEParameters();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cuda/tests/TestCudaGBSAOBCForce.cpp
View file @ 2f553a66
......@@ -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-2014 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,245 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of GBSAOBCForce.
*/
#include "CudaTests.h"
#include "TestGBSAOBCForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include "openmm/NonbondedForce.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CudaPlatform platform;
const double TOL = 1e-5;
void testSingleParticle() {
System system;
system.addParticle(2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* gbsa = new GBSAOBCForce();
NonbondedForce* nonbonded = new NonbondedForce();
gbsa->addParticle(0.5, 0.15, 1);
nonbonded->addParticle(0.5, 1, 0);
system.addForce(gbsa);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double bornRadius = 0.15-0.009; // dielectric offset
double eps0 = EPSILON0;
double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
double extendedRadius = 0.15+0.14; // probe radius
double nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
// Change the parameters and see if it is still correct.
gbsa->setParticleParameters(0, 0.4, 0.25, 1);
gbsa->updateParametersInContext(context);
state = context.getState(State::Energy);
bornRadius = 0.25-0.009; // dielectric offset
bornEnergy = (-0.4*0.4/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
extendedRadius = 0.25+0.14;
nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.25/bornRadius, 6.0);
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
}
void testGlobalSettings() {
System system;
system.addParticle(2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* forceField = new GBSAOBCForce();
forceField->addParticle(0.5, 0.15, 1);
const double soluteDielectric = 2.1;
const double solventDielectric = 35.0;
const double surfaceAreaEnergy = 0.75;
forceField->setSoluteDielectric(soluteDielectric);
forceField->setSolventDielectric(solventDielectric);
forceField->setSurfaceAreaEnergy(surfaceAreaEnergy);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double bornRadius = 0.15-0.009; // dielectric offset
double eps0 = EPSILON0;
double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/soluteDielectric-1.0/solventDielectric)/bornRadius;
double extendedRadius = 0.15+0.14; // probe radius
double nonpolarEnergy = 4*PI_M*surfaceAreaEnergy*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
}
void testCutoffAndPeriodic() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* gbsa = new GBSAOBCForce();
NonbondedForce* nonbonded = new NonbondedForce();
gbsa->addParticle(-1, 0.15, 1);
nonbonded->addParticle(-1, 1, 0);
gbsa->addParticle(1, 0.15, 1);
nonbonded->addParticle(1, 1, 0);
const double cutoffDistance = 3.0;
const double boxSize = 10.0;
nonbonded->setCutoffDistance(cutoffDistance);
gbsa->setCutoffDistance(cutoffDistance);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(gbsa);
system.addForce(nonbonded);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
// Calculate the forces for both cutoff and periodic with two different atom positions.
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
Context context(system, integrator, platform);
context.setPositions(positions);
State state1 = context.getState(State::Forces);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
context.reinitialize();
context.setPositions(positions);
State state2 = context.getState(State::Forces);
positions[1][0]+= boxSize;
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
context.reinitialize();
context.setPositions(positions);
State state3 = context.getState(State::Forces);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
context.reinitialize();
context.setPositions(positions);
State state4 = context.getState(State::Forces);
// All forces should be identical, exception state3 which should be zero.
ASSERT_EQUAL_VEC(state1.getForces()[0], state2.getForces()[0], 0.01);
ASSERT_EQUAL_VEC(state1.getForces()[1], state2.getForces()[1], 0.01);
ASSERT_EQUAL_VEC(state1.getForces()[0], state4.getForces()[0], 0.01);
ASSERT_EQUAL_VEC(state1.getForces()[1], state4.getForces()[1], 0.01);
ASSERT_EQUAL_VEC(state3.getForces()[0], Vec3(0, 0, 0), 0.01);
ASSERT_EQUAL_VEC(state3.getForces()[1], Vec3(0, 0, 0), 0.01);
void runPlatformTests() {
}
void testForce(int numParticles, NonbondedForce::NonbondedMethod method, GBSAOBCForce::NonbondedMethod method2) {
ReferencePlatform reference;
System system;
GBSAOBCForce* gbsa = new GBSAOBCForce();
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
double charge = i%2 == 0 ? -1 : 1;
gbsa->addParticle(charge, 0.15, 1);
nonbonded->addParticle(charge, 1, 0);
}
nonbonded->setNonbondedMethod(method);
gbsa->setNonbondedMethod(method2);
nonbonded->setCutoffDistance(3.0);
gbsa->setCutoffDistance(3.0);
int grid = (int) floor(0.5+pow(numParticles, 1.0/3.0));
if (method == NonbondedForce::CutoffPeriodic) {
double boxSize = (grid+1)*1.1;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
}
system.addForce(gbsa);
system.addForce(nonbonded);
LangevinIntegrator integrator1(0, 0.1, 0.01);
LangevinIntegrator integrator2(0, 0.1, 0.01);
Context context(system, integrator1, platform);
Context refContext(system, integrator2, reference);
// Set random (but uniformly distributed) positions for all the particles.
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < grid; i++)
for (int j = 0; j < grid; j++)
for (int k = 0; k < grid; k++)
positions[i*grid*grid+j*grid+k] = Vec3(i*1.1, j*1.1, k*1.1);
for (int i = 0; i < numParticles; ++i)
positions[i] = positions[i] + Vec3(0.5*genrand_real2(sfmt), 0.5*genrand_real2(sfmt), 0.5*genrand_real2(sfmt));
context.setPositions(positions);
refContext.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
State refState = refContext.getState(State::Forces | State::Energy);
// Make sure the CUDA and Reference platforms agree.
double norm = 0.0;
double diff = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = state.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
Vec3 delta = f-refState.getForces()[i];
diff += delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
}
norm = std::sqrt(norm);
diff = std::sqrt(diff);
ASSERT_EQUAL_TOL(0.0, diff, 0.001*norm);
ASSERT_EQUAL_TOL(state.getPotentialEnergy(), refState.getPotentialEnergy(), 1e-3);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
// (This doesn't work with cutoffs, since the energy changes discontinuously at the cutoff distance.)
if (method == NonbondedForce::NoCutoff)
{
const double delta = 0.3;
double step = 0.5*delta/norm;
vector<Vec3> positions2(numParticles), positions3(numParticles);
for (int i = 0; i < numParticles; ++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);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/delta, 1e-2)
}
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testSingleParticle();
testGlobalSettings();
testCutoffAndPeriodic();
for (int i = 5; i < 11; i++) {
testForce(i*i*i, NonbondedForce::NoCutoff, GBSAOBCForce::NoCutoff);
testForce(i*i*i, NonbondedForce::CutoffNonPeriodic, GBSAOBCForce::CutoffNonPeriodic);
testForce(i*i*i, NonbondedForce::CutoffPeriodic, GBSAOBCForce::CutoffPeriodic);
}
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/cuda/tests/TestCudaLangevinIntegrator.cpp
View file @ 2f553a66
......@@ -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,256 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of LangevinIntegrator.
*/
#include "CudaTests.h"
#include "TestLangevinIntegrator.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/LangevinIntegrator.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);
LangevinIntegrator integrator(0, 0.1, 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 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);
}
// Not set the friction to a tiny value and see if it conserves energy.
integrator.setFriction(5e-5);
context.setPositions(positions);
State state = context.getState(State::Energy);
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.01);
integrator.step(1);
}
}
void testTemperature() {
const int numParticles = 8;
const double temp = 100.0;
System system;
LangevinIntegrator integrator(temp, 2.0, 0.01);
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);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 10000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
}
ke /= 10000;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 6/std::sqrt(10000.0));
}
void testConstraints() {
const int numParticles = 8;
const int numConstraints = 5;
const double temp = 100.0;
System system;
LangevinIntegrator integrator(temp, 2.0, 0.01);
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.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions);
for (int j = 0; j < numConstraints; ++j) {
int particle1, particle2;
double distance;
system.getConstraintParameters(j, particle1, particle2, distance);
Vec3 p1 = state.getPositions()[particle1];
Vec3 p2 = state.getPositions()[particle2];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
ASSERT_EQUAL_TOL(distance, dist, 1e-4);
}
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);
LangevinIntegrator 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;
LangevinIntegrator integrator(temp, 2.0, 0.01);
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]);
}
}
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testSingleBond();
testTemperature();
testConstraints();
testConstrainedMasslessParticles();
testRandomSeed();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cuda/tests/TestCudaLocalEnergyMinimizer.cpp
View file @ 2f553a66
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
......@@ -7,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) 2010-2015 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -30,187 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/AssertionUtilities.h"
#include "CudaPlatform.h"
#include "openmm/Context.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/LocalEnergyMinimizer.h"
#include "openmm/NonbondedForce.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;
void testHarmonicBonds() {
const int numParticles = 10;
System system;
HarmonicBondForce* bonds = new HarmonicBondForce();
system.addForce(bonds);
// Create a chain of particles connected by harmonic bonds.
vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
positions[i] = Vec3(i, 0, 0);
if (i > 0)
bonds->addBond(i-1, i, 1+0.1*i, 1);
}
// Minimize it and check that all bonds are at their equilibrium distances.
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
LocalEnergyMinimizer::minimize(context, 1e-5);
State state = context.getState(State::Positions);
for (int i = 1; i < numParticles; i++) {
Vec3 delta = state.getPositions()[i]-state.getPositions()[i-1];
ASSERT_EQUAL_TOL(1+0.1*i, sqrt(delta.dot(delta)), 1e-4);
}
}
void testLargeSystem() {
const int numMolecules = 25;
const int numParticles = numMolecules*2;
const double cutoff = 2.0;
const double boxSize = 4.0;
const double tolerance = 10;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
system.addForce(nonbonded);
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
nonbonded->addParticle(-1.0, 0.2, 0.2);
nonbonded->addParticle(1.0, 0.2, 0.2);
positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
system.addConstraint(2*i, 2*i+1, 1.0);
}
// Minimize it and verify that the energy has decreased.
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
State initialState = context.getState(State::Forces | State::Energy);
LocalEnergyMinimizer::minimize(context, tolerance);
State finalState = context.getState(State::Forces | State::Energy | State::Positions);
ASSERT(finalState.getPotentialEnergy() < initialState.getPotentialEnergy());
// Compute the force magnitude, subtracting off any component parallel to a constraint, and
// check that it satisfies the requested tolerance.
double forceNorm = 0.0;
for (int i = 0; i < numParticles; i += 2) {
Vec3 dir = finalState.getPositions()[i+1]-finalState.getPositions()[i];
double distance = sqrt(dir.dot(dir));
dir *= 1.0/distance;
Vec3 f = finalState.getForces()[i];
f -= dir*dir.dot(f);
forceNorm += f.dot(f);
f = finalState.getForces()[i+1];
f -= dir*dir.dot(f);
forceNorm += f.dot(f);
}
forceNorm = sqrt(forceNorm/(5*numMolecules));
ASSERT(forceNorm < 2*tolerance);
}
void testVirtualSites() {
const int numMolecules = 25;
const int numParticles = numMolecules*3;
const double cutoff = 2.0;
const double boxSize = 4.0;
const double tolerance = 10;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setCutoffDistance(cutoff);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
system.addForce(nonbonded);
// Create a cloud of molecules.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(numParticles);
for (int i = 0; i < numMolecules; i++) {
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
nonbonded->addParticle(-1.0, 0.2, 0.2);
nonbonded->addParticle(0.5, 0.2, 0.2);
nonbonded->addParticle(0.5, 0.2, 0.2);
positions[3*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions[3*i+1] = Vec3(positions[3*i][0]+1.0, positions[3*i][1], positions[3*i][2]);
positions[3*i+2] = Vec3();
system.addConstraint(3*i, 3*i+1, 1.0);
system.setVirtualSite(3*i+2, new TwoParticleAverageSite(3*i, 3*i+1, 0.5, 0.5));
}
// Minimize it and verify that the energy has decreased.
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
context.applyConstraints(1e-5);
State initialState = context.getState(State::Forces | State::Energy);
LocalEnergyMinimizer::minimize(context, tolerance);
State finalState = context.getState(State::Forces | State::Energy | State::Positions);
ASSERT(finalState.getPotentialEnergy() < initialState.getPotentialEnergy());
// Compute the force magnitude, subtracting off any component parallel to a constraint, and
// check that it satisfies the requested tolerance.
double forceNorm = 0.0;
for (int i = 0; i < numParticles; i += 3) {
Vec3 dir = finalState.getPositions()[i+1]-finalState.getPositions()[i];
double distance = sqrt(dir.dot(dir));
dir *= 1.0/distance;
Vec3 f = finalState.getForces()[i];
f -= dir*dir.dot(f);
forceNorm += f.dot(f);
f = finalState.getForces()[i+1];
f -= dir*dir.dot(f);
forceNorm += f.dot(f);
// Check the virtual site location.
ASSERT_EQUAL_VEC((finalState.getPositions()[i+1]+finalState.getPositions()[i])*0.5, finalState.getPositions()[i+2], 1e-5);
}
forceNorm = sqrt(forceNorm/(5*numMolecules));
ASSERT(forceNorm < 2*tolerance);
}
#include "CudaTests.h"
#include "TestLocalEnergyMinimizer.h"
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testHarmonicBonds();
testLargeSystem();
testVirtualSites();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cuda/tests/TestCudaMonteCarloAnisotropicBarostat.cpp
View file @ 2f553a66
......@@ -6,8 +6,8 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2013 Stanford University and the Authors. *
* Authors: Peter Eastman, Lee-Ping Wang *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
......@@ -29,449 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of MonteCarloAnisotropicBarostat.
*/
#include "CudaTests.h"
#include "TestMonteCarloAnisotropicBarostat.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/MonteCarloBarostat.h"
#include "openmm/MonteCarloAnisotropicBarostat.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include "SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CudaPlatform platform;
void testIdealGas() {
const int numParticles = 64;
const int frequency = 10;
const int steps = 1000;
const double pressure = 1.5;
const double pressureInMD = pressure*(AVOGADRO*1e-25);
const double temp[] = {300.0, 600.0, 1000.0};
const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
const double initialLength = std::pow(initialVolume, 1.0/3.0);
// Create a gas of noninteracting particles.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
}
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], true, true, true, frequency);
system.addForce(barostat);
// Test it for three different temperatures.
for (int i = 0; i < 3; i++) {
barostat->setTemperature(temp[i]);
LangevinIntegrator integrator(temp[i], 0.1, 0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the volume is correct.
double volume = 0.0;
for (int j = 0; j < steps; ++j) {
Vec3 box[3];
context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
volume += box[0][0]*box[1][1]*box[2][2];
integrator.step(frequency);
}
volume /= steps;
double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
}
}
void testIdealGasAxis(int axis) {
// Test scaling just one axis.
const int numParticles = 64;
const int frequency = 10;
const int steps = 1000;
const double pressure = 1.5;
const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
const double temp[] = {300.0, 600.0, 1000.0};
const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
const double initialLength = std::pow(initialVolume, 1.0/3.0);
const bool scaleX = (axis == 0);
const bool scaleY = (axis == 1);
const bool scaleZ = (axis == 2);
double boxX;
double boxY;
double boxZ;
// Create a gas of noninteracting particles.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
}
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], scaleX, scaleY, scaleZ, frequency);
system.addForce(barostat);
// Test it for three different temperatures.
for (int i = 0; i < 3; i++) {
barostat->setTemperature(temp[i]);
LangevinIntegrator integrator(temp[i], 0.1, 0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the volume is correct.
double volume = 0.0;
for (int j = 0; j < steps; ++j) {
Vec3 box[3];
context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
boxX = box[0][0];
boxY = box[1][1];
boxZ = box[2][2];
volume += box[0][0]*box[1][1]*box[2][2];
integrator.step(frequency);
}
volume /= steps;
double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
if (!scaleX) {
ASSERT(boxX == initialLength);
}
if (!scaleY) {
ASSERT(boxY == 0.5*initialLength);
}
if (!scaleZ) {
ASSERT(boxZ == 2*initialLength);
}
}
}
void testRandomSeed() {
const int numParticles = 8;
const double temp = 100.0;
const double pressure = 1.5;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(8, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 8));
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
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);
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp, true, true, true, 1);
system.addForce(barostat);
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.
barostat->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.
barostat->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 testTriclinic() {
const int numParticles = 64;
const int frequency = 10;
const int steps = 1000;
const double pressure = 1.5;
const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
const double temperature = 300.0;
const double initialVolume = numParticles*BOLTZ*temperature/pressureInMD;
const double initialLength = std::pow(initialVolume, 1.0/3.0);
// Create a gas of noninteracting particles.
System system;
Vec3 initialBox[3];
initialBox[0] = Vec3(initialLength, 0, 0);
initialBox[1] = Vec3(0.2*initialLength, initialLength, 0);
initialBox[2] = Vec3(0.1*initialLength, 0.3*initialLength, initialLength);
system.setDefaultPeriodicBoxVectors(initialBox[0], initialBox[1], initialBox[2]);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
positions[i] = Vec3(initialLength*genrand_real2(sfmt), initialLength*genrand_real2(sfmt), initialLength*genrand_real2(sfmt));
}
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temperature, true, true, true, frequency);
system.addForce(barostat);
// Run a simulation
LangevinIntegrator integrator(temperature, 0.1, 0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the volume is correct.
double volume = 0.0;
for (int j = 0; j < steps; ++j) {
Vec3 box[3];
context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
volume += box[0][0]*box[1][1]*box[2][2];
integrator.step(frequency);
}
volume /= steps;
double expected = (numParticles+1)*BOLTZ*temperature/pressureInMD;
ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
// Make sure the box vectors have been scaled consistently.
State state = context.getState(State::Positions);
Vec3 box[3];
state.getPeriodicBoxVectors(box[0], box[1], box[2]);
double xscale = box[2][0]/(0.1*initialLength);
double yscale = box[2][1]/(0.3*initialLength);
double zscale = box[2][2]/(1.0*initialLength);
for (int i = 0; i < 3; i++) {
ASSERT_EQUAL_VEC(Vec3(xscale*initialBox[i][0], yscale*initialBox[i][1], zscale*initialBox[i][2]), box[i], 1e-5);
}
// The barostat should have put all particles inside the first periodic box. One integration step
// has happened since then, so they may have moved slightly outside it.
for (int i = 0; i < numParticles; i++) {
Vec3 pos = state.getPositions()[i];
ASSERT(pos[2]/box[2][2] > -1 && pos[2]/box[2][2] < 2);
pos -= box[2]*floor(pos[2]/box[2][2]);
ASSERT(pos[1]/box[1][1] > -1 && pos[1]/box[1][1] < 2);
pos -= box[1]*floor(pos[1]/box[1][1]);
ASSERT(pos[0]/box[0][0] > -1 && pos[0]/box[0][0] < 2);
}
}
/**
* Run a constant pressure simulation on an anisotropic Einstein crystal
* using isotropic and anisotropic barostats. There are a total of 15 simulations:
*
* 1) 3 pressures: 9.0, 10.0, 11.0 bar, for each of the following groups:
* 2) 3 groups of simulations that scale just one axis: x, y, z
* 3) 1 group of simulations that scales all three axes in the anisotropic barostat
* 4) 1 group of simulations that scales all three axes in the isotropic barostat
*
* Results that we will check:
*
* a) In each group of simulations, the volume should decrease with increasing pressure
* b) In the three simulation groups that scale just one axis, the compressibility (i.e. incremental volume change
* with increasing pressure) should go like kx > ky > kz (because the spring constant is largest in the z-direction)
* c) The anisotropic barostat should produce the same result as the isotropic barostat when all three axes are scaled
*/
void testEinsteinCrystal() {
const int numParticles = 64;
const int frequency = 2;
const int equil = 10000;
const int steps = 5000;
const double pressure = 10.0;
const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
const double temp = 300.0; // Only test one temperature since we're looking at three pressures.
const double pres3[] = {2.0, 8.0, 15.0};
const double initialVolume = numParticles*BOLTZ*temp/pressureInMD;
const double initialLength = std::pow(initialVolume, 1.0/3.0);
vector<double> initialPositions(3);
vector<double> results;
// Run four groups of anisotropic simulations; scaling just x, y, z, then all three.
for (int a = 0; a < 4; a++) {
// Test barostat for three different pressures.
for (int p = 0; p < 3; p++) {
// Create a system of noninteracting particles attached by harmonic springs to their initial positions.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, initialLength, 0), Vec3(0, 0, initialLength));
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
// Anisotropic force constants.
CustomExternalForce* force = new CustomExternalForce("0.005*(x-x0)^2 + 0.01*(y-y0)^2 + 0.02*(z-z0)^2");
force->addPerParticleParameter("x0");
force->addPerParticleParameter("y0");
force->addPerParticleParameter("z0");
NonbondedForce* nb = new NonbondedForce();
nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
positions[i] = Vec3(((i/16)%4+0.5)*initialLength/4, ((i/4)%4+0.5)*initialLength/4, (i%4+0.5)*initialLength/4);
initialPositions[0] = positions[i][0];
initialPositions[1] = positions[i][1];
initialPositions[2] = positions[i][2];
force->addParticle(i, initialPositions);
nb->addParticle(0, initialLength/6, 0.1);
}
system.addForce(force);
system.addForce(nb);
// Create the barostat.
MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pres3[p], pres3[p], pres3[p]), temp, (a==0||a==3), (a==1||a==3), (a==2||a==3), frequency);
system.addForce(barostat);
barostat->setTemperature(temp);
LangevinIntegrator integrator(temp, 0.1, 0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
// Let it equilibrate.
integrator.step(equil);
// Now run it for a while and see if the volume is correct.
double volume = 0.0;
for (int j = 0; j < steps; ++j) {
Vec3 box[3];
context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
volume += box[0][0]*box[1][1]*box[2][2];
integrator.step(frequency);
}
volume /= steps;
results.push_back(volume);
}
}
for (int p = 0; p < 3; p++) {
// Create a system of noninteracting particles attached by harmonic springs to their initial positions.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, initialLength, 0), Vec3(0, 0, initialLength));
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
// Anisotropic force constants.
CustomExternalForce* force = new CustomExternalForce("0.005*(x-x0)^2 + 0.01*(y-y0)^2 + 0.02*(z-z0)^2");
force->addPerParticleParameter("x0");
force->addPerParticleParameter("y0");
force->addPerParticleParameter("z0");
NonbondedForce* nb = new NonbondedForce();
nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
positions[i] = Vec3(((i/16)%4+0.5)*initialLength/4, ((i/4)%4+0.5)*initialLength/4, (i%4+0.5)*initialLength/4);
initialPositions[0] = positions[i][0];
initialPositions[1] = positions[i][1];
initialPositions[2] = positions[i][2];
force->addParticle(i, initialPositions);
nb->addParticle(0, initialLength/6, 0.1);
}
system.addForce(force);
system.addForce(nb);
// Create the barostat.
MonteCarloBarostat* barostat = new MonteCarloBarostat(pres3[p], temp, frequency);
system.addForce(barostat);
barostat->setTemperature(temp);
LangevinIntegrator integrator(temp, 0.1, 0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
// Let it equilibrate.
integrator.step(equil);
// Now run it for a while and see if the volume is correct.
double volume = 0.0;
for (int j = 0; j < steps; ++j) {
Vec3 box[3];
context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
volume += box[0][0]*box[1][1]*box[2][2];
integrator.step(frequency);
}
volume /= steps;
results.push_back(volume);
}
// Check to see if volumes decrease with increasing pressure.
ASSERT_USUALLY_TRUE(results[0] > results[1]);
ASSERT_USUALLY_TRUE(results[1] > results[2]);
ASSERT_USUALLY_TRUE(results[3] > results[4]);
ASSERT_USUALLY_TRUE(results[4] > results[5]);
ASSERT_USUALLY_TRUE(results[6] > results[7]);
ASSERT_USUALLY_TRUE(results[7] > results[8]);
// Check to see if incremental volume changes with increasing pressure go like kx > ky > kz.
ASSERT_USUALLY_TRUE((results[0] - results[1]) > (results[3] - results[4]));
ASSERT_USUALLY_TRUE((results[1] - results[2]) > (results[4] - results[5]));
ASSERT_USUALLY_TRUE((results[3] - results[4]) > (results[6] - results[7]));
ASSERT_USUALLY_TRUE((results[4] - results[5]) > (results[7] - results[8]));
// Check to see if the volumes are equal for isotropic and anisotropic (all axis).
ASSERT_USUALLY_EQUAL_TOL(results[9], results[12], 3/std::sqrt((double) steps));
ASSERT_USUALLY_EQUAL_TOL(results[10], results[13], 3/std::sqrt((double) steps));
ASSERT_USUALLY_EQUAL_TOL(results[11], results[14], 3/std::sqrt((double) steps));
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testIdealGas();
testIdealGasAxis(0);
testIdealGasAxis(1);
testIdealGasAxis(2);
testRandomSeed();
testTriclinic();
//testEinsteinCrystal();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cuda/tests/TestCudaMonteCarloBarostat.cpp
View file @ 2f553a66
......@@ -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,261 +29,9 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of MonteCarloBarostat.
*/
#include "CudaTests.h"
#include "TestMonteCarloBarostat.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/MonteCarloBarostat.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include "SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CudaPlatform platform;
void testChangingBoxSize() {
System system;
system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 5, 0), Vec3(0, 0, 6));
system.addParticle(1.0);
NonbondedForce* nb = new NonbondedForce();
nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nb->setCutoffDistance(2.0);
nb->addParticle(1, 0.5, 0.5);
system.addForce(nb);
LangevinIntegrator integrator(300.0, 1.0, 0.01);
Context context(system, integrator, platform);
vector<Vec3> positions;
positions.push_back(Vec3());
context.setPositions(positions);
Vec3 x, y, z;
context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
ASSERT_EQUAL_VEC(Vec3(4, 0, 0), x, 0);
ASSERT_EQUAL_VEC(Vec3(0, 5, 0), y, 0);
ASSERT_EQUAL_VEC(Vec3(0, 0, 6), z, 0);
context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 9));
context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
ASSERT_EQUAL_VEC(Vec3(7, 0, 0), x, 0);
ASSERT_EQUAL_VEC(Vec3(0, 8, 0), y, 0);
ASSERT_EQUAL_VEC(Vec3(0, 0, 9), z, 0);
// Shrinking the box too small should produce an exception.
context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 3.9, 0), Vec3(0, 0, 9));
bool ok = true;
try {
context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
ok = false;
}
catch (exception& ex) {
}
ASSERT(ok);
}
void testIdealGas() {
const int numParticles = 64;
const int frequency = 10;
const int steps = 1000;
const double pressure = 1.5;
const double pressureInMD = pressure*(AVOGADRO*1e-25);
const double temp[] = {300.0, 600.0, 1000.0};
const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
const double initialLength = std::pow(initialVolume, 1.0/3.0);
// Create a gas of noninteracting particles.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
}
MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency);
system.addForce(barostat);
// Test it for three different temperatures.
for (int i = 0; i < 3; i++) {
barostat->setTemperature(temp[i]);
LangevinIntegrator integrator(temp[i], 0.1, 0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the volume is correct.
double volume = 0.0;
for (int j = 0; j < steps; ++j) {
Vec3 box[3];
context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
volume += box[0][0]*box[1][1]*box[2][2];
ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
integrator.step(frequency);
}
volume /= steps;
double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
}
}
void testRandomSeed() {
const int numParticles = 8;
const double temp = 100.0;
const double pressure = 1.5;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(8, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 8));
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
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);
MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, 1);
system.addForce(barostat);
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.
barostat->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.
barostat->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 testWater() {
const int gridSize = 8;
const int numMolecules = gridSize*gridSize*gridSize;
const int frequency = 10;
const int steps = 400;
const double temp = 273.15;
const double pressure = 3;
const double spacing = 0.32;
const double angle = 109.47*M_PI/180;
const double dOH = 0.1;
const double dHH = dOH*2*std::sin(0.5*angle);
// Create a box of SPC water molecules.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(gridSize*spacing, 0, 0), Vec3(0, gridSize*spacing, 0), Vec3(0, 0, gridSize*spacing));
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setUseDispersionCorrection(true);
vector<Vec3> positions;
Vec3 offset1(dOH, 0, 0);
Vec3 offset2(dOH*std::cos(angle), dOH*std::sin(angle), 0);
for (int i = 0; i < gridSize; ++i) {
for (int j = 0; j < gridSize; ++j) {
for (int k = 0; k < gridSize; ++k) {
int firstParticle = system.getNumParticles();
system.addParticle(16.0);
system.addParticle(1.0);
system.addParticle(1.0);
nonbonded->addParticle(-0.82, 0.316557, 0.650194);
nonbonded->addParticle(0.41, 1, 0);
nonbonded->addParticle(0.41, 1, 0);
Vec3 pos = Vec3(spacing*i, spacing*j, spacing*k);
positions.push_back(pos);
positions.push_back(pos+offset1);
positions.push_back(pos+offset2);
system.addConstraint(firstParticle, firstParticle+1, dOH);
system.addConstraint(firstParticle, firstParticle+2, dOH);
system.addConstraint(firstParticle+1, firstParticle+2, dHH);
nonbonded->addException(firstParticle, firstParticle+1, 0, 1, 0);
nonbonded->addException(firstParticle, firstParticle+2, 0, 1, 0);
nonbonded->addException(firstParticle+1, firstParticle+2, 0, 1, 0);
}
}
}
system.addForce(nonbonded);
MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency);
system.addForce(barostat);
// Simulate it and see if the density matches the expected value (1 g/mL).
LangevinIntegrator integrator(temp, 1.0, 0.002);
Context context(system, integrator, platform);
context.setPositions(positions);
integrator.step(2000);
double volume = 0.0;
for (int j = 0; j < steps; ++j) {
Vec3 box[3];
context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
volume += box[0][0]*box[1][1]*box[2][2];
integrator.step(frequency);
}
volume /= steps;
double density = numMolecules*18/(AVOGADRO*volume*1e-21);
ASSERT_USUALLY_EQUAL_TOL(1.0, density, 0.02);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testChangingBoxSize();
testIdealGas();
testRandomSeed();
testWater();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
testWater();
}
platforms/cuda/tests/TestCudaNonbondedForce.cpp
View file @ 2f553a66
......@@ -29,785 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests all the different force terms in the CUDA implementation of NonbondedForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/internal/ContextImpl.h"
#include "CudaArray.h"
#include "CudaNonbondedUtilities.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 testCoulomb() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->addParticle(0.5, 1, 0);
forceField->addParticle(-1.5, 1, 0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double force = ONE_4PI_EPS0*(-0.75)/4.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(-0.75)/2.0, state.getPotentialEnergy(), TOL);
}
void testLJ() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->addParticle(0, 1.2, 1);
forceField->addParticle(0, 1.4, 2);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double x = 1.3/2.0;
double eps = SQRT_TWO;
double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/2.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0)), state.getPotentialEnergy(), TOL);
}
void testExclusionsAnd14() {
System system;
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < 5; ++i) {
system.addParticle(1.0);
nonbonded->addParticle(0, 1.5, 0);
}
vector<pair<int, int> > bonds;
bonds.push_back(pair<int, int>(0, 1));
bonds.push_back(pair<int, int>(1, 2));
bonds.push_back(pair<int, int>(2, 3));
bonds.push_back(pair<int, int>(3, 4));
nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
int first14, second14;
for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon;
nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
first14 = i;
if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
second14 = i;
}
system.addForce(nonbonded);
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
for (int i = 1; i < 5; ++i) {
// Test LJ forces
vector<Vec3> positions(5);
const double r = 1.0;
for (int j = 0; j < 5; ++j) {
nonbonded->setParticleParameters(j, 0, 1.5, 0);
positions[j] = Vec3(0, j, 0);
}
nonbonded->setParticleParameters(0, 0, 1.5, 1);
nonbonded->setParticleParameters(i, 0, 1.5, 1);
nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
positions[i] = Vec3(r, 0, 0);
context.reinitialize();
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double x = 1.5/r;
double eps = 1.0;
double force = 4.0*eps*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
double energy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
if (i == 3) {
force *= 0.5;
energy *= 0.5;
}
if (i < 3) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
// Test Coulomb forces
nonbonded->setParticleParameters(0, 2, 1.5, 0);
nonbonded->setParticleParameters(i, 2, 1.5, 0);
nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces();
force = ONE_4PI_EPS0*4/(r*r);
energy = ONE_4PI_EPS0*4/r;
if (i == 3) {
force /= 1.2;
energy /= 1.2;
}
if (i < 3) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
}
void testCutoff() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce();
forceField->addParticle(1.0, 1, 0);
forceField->addParticle(1.0, 1, 0);
forceField->addParticle(1.0, 1, 0);
forceField->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
const double cutoff = 2.9;
forceField->setCutoffDistance(cutoff);
const double eps = 50.0;
forceField->setReactionFieldDielectric(eps);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(0, 2, 0);
positions[2] = Vec3(0, 3, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force1 = ONE_4PI_EPS0*(1.0)*(0.25-2.0*krf*2.0);
const double force2 = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL);
const double energy1 = ONE_4PI_EPS0*(1.0)*(0.5+krf*4.0-crf);
const double energy2 = ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf);
ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL);
}
void testCutoff14() {
System system;
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
for (int i = 0; i < 5; ++i) {
system.addParticle(1.0);
nonbonded->addParticle(0, 1.5, 0);
}
const double cutoff = 3.5;
nonbonded->setCutoffDistance(cutoff);
const double eps = 30.0;
nonbonded->setReactionFieldDielectric(eps);
vector<pair<int, int> > bonds;
bonds.push_back(pair<int, int>(0, 1));
bonds.push_back(pair<int, int>(1, 2));
bonds.push_back(pair<int, int>(2, 3));
bonds.push_back(pair<int, int>(3, 4));
nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
int first14, second14;
for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon;
nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
first14 = i;
if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
second14 = i;
}
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(5);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(2, 0, 0);
positions[3] = Vec3(3, 0, 0);
positions[4] = Vec3(4, 0, 0);
for (int i = 1; i < 5; ++i) {
// Test LJ forces
nonbonded->setParticleParameters(0, 0, 1.5, 1);
for (int j = 1; j < 5; ++j)
nonbonded->setParticleParameters(j, 0, 1.5, 0);
nonbonded->setParticleParameters(i, 0, 1.5, 1);
nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
context.reinitialize();
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double r = positions[i][0];
double x = 1.5/r;
double e = 1.0;
double force = 4.0*e*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
double energy = 4.0*e*(std::pow(x, 12.0)-std::pow(x, 6.0));
if (i == 3) {
force *= 0.5;
energy *= 0.5;
}
if (i < 3 || r > cutoff) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
// Test Coulomb forces
const double q = 0.7;
nonbonded->setParticleParameters(0, q, 1.5, 0);
nonbonded->setParticleParameters(i, q, 1.5, 0);
nonbonded->setExceptionParameters(first14, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces();
force = ONE_4PI_EPS0*q*q/(r*r);
energy = ONE_4PI_EPS0*q*q/r;
if (i == 3) {
force /= 1.2;
energy /= 1.2;
}
if (i < 3 || r > cutoff) {
force = 0;
energy = 0;
}
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
}
}
void testPeriodic() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addException(0, 1, 0.0, 1.0, 0.0);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
positions[2] = Vec3(3, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
const double eps = 78.3;
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
}
void testTriclinic() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
Vec3 a(3.1, 0, 0);
Vec3 b(0.4, 3.5, 0);
Vec3 c(-0.1, -0.5, 4.0);
system.setDefaultPeriodicBoxVectors(a, b, c);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 1.5;
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const double eps = 78.3;
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
for (int iteration = 0; iteration < 50; iteration++) {
// Generate random positions for the two particles.
positions[0] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
positions[1] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
context.setPositions(positions);
// Loop over all possible periodic copies and find the nearest one.
Vec3 delta;
double distance2 = 100.0;
for (int i = -1; i < 2; i++)
for (int j = -1; j < 2; j++)
for (int k = -1; k < 2; k++) {
Vec3 d = positions[1]-positions[0]+a*i+b*j+c*k;
if (d.dot(d) < distance2) {
delta = d;
distance2 = d.dot(d);
}
}
double distance = sqrt(distance2);
// See if the force and energy are correct.
State state = context.getState(State::Forces | State::Energy);
if (distance >= cutoff) {
ASSERT_EQUAL(0.0, state.getPotentialEnergy());
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getForces()[0], 0);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getForces()[1], 0);
}
else {
const Vec3 force = delta*ONE_4PI_EPS0*(-1.0/(distance*distance*distance)+2.0*krf);
ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(1.0/distance+krf*distance*distance-crf), state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_VEC(force, state.getForces()[0], TOL);
ASSERT_EQUAL_VEC(-force, state.getForces()[1], TOL);
}
}
}
void testLargeSystem() {
const int numMolecules = 600;
const int numParticles = numMolecules*2;
const double cutoff = 2.0;
const double boxSize = 20.0;
const double tol = 2e-3;
ReferencePlatform reference;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
NonbondedForce* nonbonded = new NonbondedForce();
HarmonicBondForce* bonds = new HarmonicBondForce();
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numMolecules; i++) {
if (i < numMolecules/2) {
nonbonded->addParticle(-1.0, 0.2, 0.1);
nonbonded->addParticle(1.0, 0.1, 0.1);
}
else {
nonbonded->addParticle(-1.0, 0.2, 0.2);
nonbonded->addParticle(1.0, 0.1, 0.2);
}
positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
velocities[2*i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
velocities[2*i+1] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
bonds->addBond(2*i, 2*i+1, 1.0, 0.1);
nonbonded->addException(2*i, 2*i+1, 0.0, 0.15, 0.0);
}
// Try with cutoffs but not periodic boundary conditions, and make sure the cl and Reference
// platforms agree.
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
system.addForce(bonds);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context cuContext(system, integrator1, platform);
Context referenceContext(system, integrator2, reference);
cuContext.setPositions(positions);
cuContext.setVelocities(velocities);
referenceContext.setPositions(positions);
referenceContext.setVelocities(velocities);
State cuState = cuContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(cuState.getPositions()[i], referenceState.getPositions()[i], tol);
ASSERT_EQUAL_VEC(cuState.getVelocities()[i], referenceState.getVelocities()[i], tol);
ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
}
ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
// Now do the same thing with periodic boundary conditions.
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
cuContext.reinitialize();
referenceContext.reinitialize();
cuContext.setPositions(positions);
cuContext.setVelocities(velocities);
referenceContext.setPositions(positions);
referenceContext.setVelocities(velocities);
cuState = cuContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
referenceState = referenceContext.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++) {
double dx = cuState.getPositions()[i][0]-referenceState.getPositions()[i][0];
double dy = cuState.getPositions()[i][1]-referenceState.getPositions()[i][1];
double dz = cuState.getPositions()[i][2]-referenceState.getPositions()[i][2];
ASSERT_EQUAL_TOL(fmod(cuState.getPositions()[i][0]-referenceState.getPositions()[i][0], boxSize), 0, tol);
ASSERT_EQUAL_TOL(fmod(cuState.getPositions()[i][1]-referenceState.getPositions()[i][1], boxSize), 0, tol);
ASSERT_EQUAL_TOL(fmod(cuState.getPositions()[i][2]-referenceState.getPositions()[i][2], boxSize), 0, tol);
ASSERT_EQUAL_VEC(cuState.getVelocities()[i], referenceState.getVelocities()[i], tol);
ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
}
ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
}
/*
void testBlockInteractions(bool periodic) {
const int blockSize = CudaContext::TileSize;
const int numBlocks = 100;
const int numParticles = blockSize*numBlocks;
const double cutoff = 1.0;
const double boxSize = (periodic ? 5.1 : 1.1);
System system;
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
nonbonded->addParticle(1.0, 0.2, 0.2);
positions[i] = Vec3(boxSize*(3*genrand_real2(sfmt)-1), boxSize*(3*genrand_real2(sfmt)-1), boxSize*(3*genrand_real2(sfmt)-1));
}
nonbonded->setNonbondedMethod(periodic ? NonbondedForce::CutoffPeriodic : NonbondedForce::CutoffNonPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
Context context(system, integrator, platform);
context.setPositions(positions);
ContextImpl* contextImpl = *reinterpret_cast<ContextImpl**>(&context);
CudaPlatform::PlatformData& data = *static_cast<CudaPlatform::PlatformData*>(contextImpl->getPlatformData());
CudaContext& cuContext = *data.contexts[0];
CudaNonbondedUtilities& nb = cuContext.getNonbondedUtilities();
State state = context.getState(State::Positions | State::Velocities | State::Forces);
nb.updateNeighborListSize();
state = context.getState(State::Positions | State::Velocities | State::Forces);
// Verify that the bounds of each block were calculated correctly.
vector<double4> posq(cuContext.getPosq().getSize());
vector<double4> blockCenters(numBlocks);
vector<double4> blockBoundingBoxes(numBlocks);
if (cuContext.getUseDoublePrecision()) {
cuContext.getPosq().download(posq);
nb.getBlockCenters().download(blockCenters);
nb.getBlockBoundingBoxes().download(blockBoundingBoxes);
}
else {
vector<float4> posqf(cuContext.getPosq().getSize());
vector<float4> blockCentersf(numBlocks);
vector<float4> blockBoundingBoxesf(numBlocks);
cuContext.getPosq().download(posqf);
nb.getBlockCenters().download(blockCentersf);
nb.getBlockBoundingBoxes().download(blockBoundingBoxesf);
for (int i = 0; i < numParticles; i++)
posq[i] = make_double4(posqf[i].x, posqf[i].y, posqf[i].z, posqf[i].w);
for (int i = 0; i < numBlocks; i++) {
blockCenters[i] = make_double4(blockCentersf[i].x, blockCentersf[i].y, blockCentersf[i].z, blockCentersf[i].w);
blockBoundingBoxes[i] = make_double4(blockBoundingBoxesf[i].x, blockBoundingBoxesf[i].y, blockBoundingBoxesf[i].z, blockBoundingBoxesf[i].w);
}
}
for (int i = 0; i < numBlocks; i++) {
double4 gridSize = blockBoundingBoxes[i];
double4 center = blockCenters[i];
if (periodic) {
ASSERT(gridSize.x < 0.5*boxSize);
ASSERT(gridSize.y < 0.5*boxSize);
ASSERT(gridSize.z < 0.5*boxSize);
}
double minx = 0.0, maxx = 0.0, miny = 0.0, maxy = 0.0, minz = 0.0, maxz = 0.0, radius = 0.0;
for (int j = 0; j < blockSize; j++) {
double4 pos = posq[i*blockSize+j];
double dx = pos.x-center.x;
double dy = pos.y-center.y;
double dz = pos.z-center.z;
if (periodic) {
dx -= floor(0.5+dx/boxSize)*boxSize;
dy -= floor(0.5+dy/boxSize)*boxSize;
dz -= floor(0.5+dz/boxSize)*boxSize;
}
ASSERT(abs(dx) < gridSize.x+TOL);
ASSERT(abs(dy) < gridSize.y+TOL);
ASSERT(abs(dz) < gridSize.z+TOL);
minx = min(minx, dx);
maxx = max(maxx, dx);
miny = min(miny, dy);
maxy = max(maxy, dy);
minz = min(minz, dz);
maxz = max(maxz, dz);
}
ASSERT_EQUAL_TOL(-minx, gridSize.x, TOL);
ASSERT_EQUAL_TOL(maxx, gridSize.x, TOL);
ASSERT_EQUAL_TOL(-miny, gridSize.y, TOL);
ASSERT_EQUAL_TOL(maxy, gridSize.y, TOL);
ASSERT_EQUAL_TOL(-minz, gridSize.z, TOL);
ASSERT_EQUAL_TOL(maxz, gridSize.z, TOL);
}
// Verify that interactions were identified correctly.
vector<unsigned int> interactionCount;
vector<ushort2> interactingTiles;
vector<unsigned int> interactionFlags;
nb.getInteractionCount().download(interactionCount);
int numWithInteractions = interactionCount[0];
vector<bool> hasInteractions(numBlocks*(numBlocks+1)/2, false);
nb.getInteractingTiles().download(interactingTiles);
nb.getInteractionFlags().download(interactionFlags);
const unsigned int dim = cuContext.getNumAtomBlocks();
for (int i = 0; i < numWithInteractions; i++) {
unsigned int x = interactingTiles[i].x;
unsigned int y = interactingTiles[i].y;
int index = (x > y ? x+y*dim-y*(y+1)/2 : y+x*dim-x*(x+1)/2);
hasInteractions[index] = true;
// Make sure this tile really should have been flagged based on bounding volumes.
double4 gridSize1 = blockBoundingBoxes[x];
double4 gridSize2 = blockBoundingBoxes[y];
double4 center1 = blockCenters[x];
double4 center2 = blockCenters[y];
double dx = center1.x-center2.x;
double dy = center1.y-center2.y;
double dz = center1.z-center2.z;
if (periodic) {
dx -= floor(0.5+dx/boxSize)*boxSize;
dy -= floor(0.5+dy/boxSize)*boxSize;
dz -= floor(0.5+dz/boxSize)*boxSize;
}
dx = max(0.0, abs(dx)-gridSize1.x-gridSize2.x);
dy = max(0.0, abs(dy)-gridSize1.y-gridSize2.y);
dz = max(0.0, abs(dz)-gridSize1.z-gridSize2.z);
ASSERT(sqrt(dx*dx+dy*dy+dz*dz) < cutoff+TOL);
// Check the interaction flags.
unsigned int flags = interactionFlags[i];
for (int atom2 = 0; atom2 < 32; atom2++) {
if ((flags & 1) == 0) {
double4 pos2 = posq[y*blockSize+atom2];
for (int atom1 = 0; atom1 < blockSize; ++atom1) {
double4 pos1 = posq[x*blockSize+atom1];
double dx = pos2.x-pos1.x;
double dy = pos2.y-pos1.y;
double dz = pos2.z-pos1.z;
if (periodic) {
dx -= floor(0.5+dx/boxSize)*boxSize;
dy -= floor(0.5+dy/boxSize)*boxSize;
dz -= floor(0.5+dz/boxSize)*boxSize;
}
ASSERT(dx*dx+dy*dy+dz*dz > cutoff*cutoff);
}
}
flags >>= 1;
}
}
// Check the tiles that did not have interactions to make sure all atoms are beyond the cutoff.
for (int i = 0; i < (int) hasInteractions.size(); i++)
if (!hasInteractions[i]) {
unsigned int y = (unsigned int) std::floor(numBlocks+0.5-std::sqrt((numBlocks+0.5)*(numBlocks+0.5)-2*i));
unsigned int x = (i-y*numBlocks+y*(y+1)/2);
if (x == y)
continue; // This block has exclusions, so it will not be in the neighbor list.
for (int atom1 = 0; atom1 < blockSize; ++atom1) {
double4 pos1 = posq[x*blockSize+atom1];
for (int atom2 = 0; atom2 < blockSize; ++atom2) {
double4 pos2 = posq[y*blockSize+atom2];
double dx = pos1.x-pos2.x;
double dy = pos1.y-pos2.y;
double dz = pos1.z-pos2.z;
if (periodic) {
dx -= floor(0.5+dx/boxSize)*boxSize;
dy -= floor(0.5+dy/boxSize)*boxSize;
dz -= floor(0.5+dz/boxSize)*boxSize;
}
ASSERT(dx*dx+dy*dy+dz*dz > cutoff*cutoff);
}
}
}
}*/
void testDispersionCorrection() {
// Create a box full of identical particles.
int gridSize = 5;
int numParticles = gridSize*gridSize*gridSize;
double boxSize = gridSize*0.7;
double cutoff = boxSize/3;
System system;
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions(numParticles);
int index = 0;
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
system.addParticle(1.0);
nonbonded->addParticle(0, 1.1, 0.5);
positions[index] = Vec3(i*boxSize/gridSize, j*boxSize/gridSize, k*boxSize/gridSize);
index++;
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
// See if the correction has the correct value.
Context context(system, integrator, platform);
context.setPositions(positions);
double energy1 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseDispersionCorrection(false);
context.reinitialize();
context.setPositions(positions);
double energy2 = context.getState(State::Energy).getPotentialEnergy();
double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
// Now modify half the particles to be different, and see if it is still correct.
int numType2 = 0;
for (int i = 0; i < numParticles; i += 2) {
nonbonded->setParticleParameters(i, 0, 1, 1);
numType2++;
}
int numType1 = numParticles-numType2;
nonbonded->updateParametersInContext(context);
energy2 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseDispersionCorrection(true);
context.reinitialize();
context.setPositions(positions);
energy1 = context.getState(State::Energy).getPotentialEnergy();
term1 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
term2 = ((numType1*(numType1+1))/2)*(0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
term1 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 12)/pow(cutoff, 9))/9;
term2 += ((numType2*(numType2+1))/2)*(1*pow(1.0, 6)/pow(cutoff, 3))/3;
double combinedSigma = 0.5*(1+1.1);
double combinedEpsilon = sqrt(1*0.5);
term1 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 12)/pow(cutoff, 9))/9;
term2 += (numType1*numType2)*(combinedEpsilon*pow(combinedSigma, 6)/pow(cutoff, 3))/3;
term1 /= (numParticles*(numParticles+1))/2;
term2 /= (numParticles*(numParticles+1))/2;
expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
}
void testChangingParameters() {
const int numMolecules = 600;
const int numParticles = numMolecules*2;
const double cutoff = 2.0;
const double boxSize = 20.0;
const double tol = 2e-3;
ReferencePlatform reference;
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numMolecules; i++) {
if (i < numMolecules/2) {
nonbonded->addParticle(-1.0, 0.2, 0.1);
nonbonded->addParticle(1.0, 0.1, 0.1);
}
else {
nonbonded->addParticle(-1.0, 0.2, 0.2);
nonbonded->addParticle(1.0, 0.1, 0.2);
}
positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
system.addConstraint(2*i, 2*i+1, 1.0);
nonbonded->addException(2*i, 2*i+1, 0.0, 0.15, 0.0);
}
nonbonded->setNonbondedMethod(NonbondedForce::PME);
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
// See if Reference and Cuda give the same forces and energies.
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context cuContext(system, integrator1, platform);
Context referenceContext(system, integrator2, reference);
cuContext.setPositions(positions);
referenceContext.setPositions(positions);
State cuState = cuContext.getState(State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
// Now modify parameters and see if they still agree.
for (int i = 0; i < numParticles; i += 5) {
double charge, sigma, epsilon;
nonbonded->getParticleParameters(i, charge, sigma, epsilon);
nonbonded->setParticleParameters(i, 1.5*charge, 1.1*sigma, 1.7*epsilon);
}
nonbonded->updateParametersInContext(cuContext);
nonbonded->updateParametersInContext(referenceContext);
cuState = cuContext.getState(State::Forces | State::Energy);
referenceState = referenceContext.getState(State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(cuState.getForces()[i], referenceState.getForces()[i], tol);
ASSERT_EQUAL_TOL(cuState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
}
#include "CudaTests.h"
#include "TestNonbondedForce.h"
void testParallelComputation(NonbondedForce::NonbondedMethod method) {
System system;
......@@ -868,61 +91,6 @@ void testParallelComputation(NonbondedForce::NonbondedMethod method) {
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
}
void testSwitchingFunction(NonbondedForce::NonbondedMethod method) {
System system;
system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(0, 1.2, 1);
nonbonded->addParticle(0, 1.4, 2);
nonbonded->setNonbondedMethod(method);
nonbonded->setCutoffDistance(2.0);
nonbonded->setUseSwitchingFunction(true);
nonbonded->setSwitchingDistance(1.5);
nonbonded->setUseDispersionCorrection(false);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
double eps = SQRT_TWO;
// Compute the interaction at various distances.
for (double r = 1.0; r < 2.5; r += 0.1) {
positions[1] = Vec3(r, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
// See if the energy is correct.
double x = 1.3/r;
double expectedEnergy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
double switchValue;
if (r <= 1.5)
switchValue = 1;
else if (r >= 2.0)
switchValue = 0;
else {
double t = (r-1.5)/0.5;
switchValue = 1+t*t*t*(-10+t*(15-t*6));
}
ASSERT_EQUAL_TOL(switchValue*expectedEnergy, state.getPotentialEnergy(), TOL);
// See if the force is the gradient of the energy.
double delta = 1e-3;
positions[1] = Vec3(r-delta, 0, 0);
context.setPositions(positions);
double e1 = context.getState(State::Energy).getPotentialEnergy();
positions[1] = Vec3(r+delta, 0, 0);
context.setPositions(positions);
double e2 = context.getState(State::Energy).getPotentialEnergy();
ASSERT_EQUAL_TOL((e2-e1)/(2*delta), state.getForces()[0][0], 1e-3);
}
}
void testReordering() {
// Check that reordering of atoms doesn't alter their positions.
......@@ -950,33 +118,9 @@ void testReordering() {
}
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testCoulomb();
testLJ();
testExclusionsAnd14();
testCutoff();
testCutoff14();
testPeriodic();
testTriclinic();
testLargeSystem();
//testBlockInteractions(false);
//testBlockInteractions(true);
testDispersionCorrection();
testChangingParameters();
testParallelComputation(NonbondedForce::NoCutoff);
testParallelComputation(NonbondedForce::Ewald);
testParallelComputation(NonbondedForce::PME);
testSwitchingFunction(NonbondedForce::CutoffNonPeriodic);
testSwitchingFunction(NonbondedForce::PME);
testReordering();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
testParallelComputation(NonbondedForce::NoCutoff);
testParallelComputation(NonbondedForce::Ewald);
testParallelComputation(NonbondedForce::PME);
testReordering();
}
platforms/cuda/tests/TestCudaSettle.cpp
View file @ 2f553a66
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
......@@ -7,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: *
* *
......@@ -30,90 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CUDA implementation of the SETTLE algorithm.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
CudaPlatform platform;
void testConstraints() {
const int numMolecules = 10;
const int numParticles = numMolecules*3;
const int numConstraints = numMolecules*3;
const double temp = 100.0;
System system;
LangevinIntegrator integrator(temp, 2.0, 0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numMolecules; ++i) {
system.addParticle(16.0);
system.addParticle(1.0);
system.addParticle(1.0);
forceField->addParticle(-0.82, 0.317, 0.65);
forceField->addParticle(0.41, 1.0, 0.0);
forceField->addParticle(0.41, 1.0, 0.0);
system.addConstraint(i*3, i*3+1, 0.1);
system.addConstraint(i*3, i*3+2, 0.1);
system.addConstraint(i*3+1, i*3+2, 0.163);
}
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 < numMolecules; ++i) {
positions[i*3] = Vec3((i%4)*0.4, (i/4)*0.4, 0);
positions[i*3+1] = positions[i*3]+Vec3(0.1, 0, 0);
positions[i*3+2] = positions[i*3]+Vec3(-0.03333, 0.09428, 0);
velocities[i*3] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
velocities[i*3+1] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
velocities[i*3+2] = 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) {
integrator.step(1);
State state = context.getState(State::Positions | 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-5);
}
}
}
#include "CudaTests.h"
#include "TestSettle.h"
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testConstraints();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
}
platforms/cuda/tests/nacl_amorph.dat deleted 100644 → 0
View file @ 4ed151e1
positions[0] = Vec3(1.066000,1.628000,0.835000);
positions[1] = Vec3(1.072000,0.428000,0.190000);
positions[2] = Vec3(0.524000,1.442000,1.160000);
positions[3] = Vec3(2.383000,1.524000,1.119000);
positions[4] = Vec3(0.390000,1.441000,0.575000);
positions[5] = Vec3(0.618000,0.399000,0.819000);
positions[6] = Vec3(1.003000,1.257000,1.543000);
positions[7] = Vec3(2.933000,1.569000,0.642000);
positions[8] = Vec3(0.849000,0.739000,0.089000);
positions[9] = Vec3(0.060000,0.794000,0.766000);
positions[10] = Vec3(1.652000,1.405000,1.010000);
positions[11] = Vec3(2.843000,1.533000,1.781000);
positions[12] = Vec3(0.952000,1.309000,0.996000);
positions[13] = Vec3(1.847000,1.402000,0.313000);
positions[14] = Vec3(2.674000,0.083000,1.691000);
positions[15] = Vec3(1.763000,2.104000,0.728000);
positions[16] = Vec3(0.914000,0.574000,0.982000);
positions[17] = Vec3(0.514000,0.078000,0.891000);
positions[18] = Vec3(0.538000,0.766000,1.110000);
positions[19] = Vec3(0.808000,0.676000,0.570000);
positions[20] = Vec3(0.178000,0.014000,0.628000);
positions[21] = Vec3(1.329000,1.333000,0.339000);
positions[22] = Vec3(1.029000,1.678000,0.503000);
positions[23] = Vec3(1.423000,1.767000,1.104000);
positions[24] = Vec3(1.966000,1.051000,0.282000);
positions[25] = Vec3(1.596000,1.971000,0.194000);
positions[26] = Vec3(1.025000,1.043000,2.809000);
positions[27] = Vec3(1.628000,2.614000,0.088000);
positions[28] = Vec3(0.440000,0.606000,0.141000);
positions[29] = Vec3(1.050000,2.821000,2.517000);
positions[30] = Vec3(0.644000,1.604000,0.770000);
positions[31] = Vec3(0.637000,0.917000,0.392000);
positions[32] = Vec3(0.611000,2.768000,0.013000);
positions[33] = Vec3(1.892000,0.660000,0.473000);
positions[34] = Vec3(1.052000,2.081000,0.982000);
positions[35] = Vec3(1.508000,2.300000,0.439000);
positions[36] = Vec3(2.617000,0.328000,1.099000);
positions[37] = Vec3(0.910000,0.040000,0.259000);
positions[38] = Vec3(1.195000,1.494000,1.202000);
positions[39] = Vec3(2.657000,0.997000,0.564000);
positions[40] = Vec3(1.465000,1.580000,0.648000);
positions[41] = Vec3(0.154000,2.538000,1.331000);
positions[42] = Vec3(0.849000,1.476000,1.365000);
positions[43] = Vec3(0.898000,0.987000,1.178000);
positions[44] = Vec3(0.958000,0.656000,1.358000);
positions[45] = Vec3(1.067000,0.934000,0.211000);
positions[46] = Vec3(1.030000,0.319000,1.281000);
positions[47] = Vec3(2.709000,0.807000,0.240000);
positions[48] = Vec3(0.837000,1.362000,0.588000);
positions[49] = Vec3(2.080000,0.791000,2.947000);
positions[50] = Vec3(0.200000,0.266000,1.474000);
positions[51] = Vec3(0.848000,0.379000,1.625000);
positions[52] = Vec3(0.637000,1.071000,0.821000);
positions[53] = Vec3(1.324000,0.757000,2.951000);
positions[54] = Vec3(2.666000,0.935000,1.373000);
positions[55] = Vec3(1.584000,1.025000,1.703000);
positions[56] = Vec3(1.699000,0.636000,0.038000);
positions[57] = Vec3(1.099000,1.644000,1.879000);
positions[58] = Vec3(2.897000,1.302000,1.522000);
positions[59] = Vec3(1.753000,0.949000,2.885000);
positions[60] = Vec3(2.502000,1.321000,0.752000);
positions[61] = Vec3(0.545000,0.193000,1.959000);
positions[62] = Vec3(1.098000,2.646000,1.706000);
positions[63] = Vec3(0.001000,1.205000,0.670000);
positions[64] = Vec3(2.997000,0.061000,1.040000);
positions[65] = Vec3(0.662000,0.828000,1.535000);
positions[66] = Vec3(1.252000,1.246000,0.780000);
positions[67] = Vec3(1.173000,0.472000,0.810000);
positions[68] = Vec3(0.124000,0.622000,2.992000);
positions[69] = Vec3(1.036000,0.883000,0.848000);
positions[70] = Vec3(1.423000,2.146000,1.340000);
positions[71] = Vec3(2.391000,1.136000,1.165000);
positions[72] = Vec3(1.189000,2.961000,0.425000);
positions[73] = Vec3(1.584000,2.500000,0.782000);
positions[74] = Vec3(0.565000,1.122000,1.240000);
positions[75] = Vec3(1.733000,1.716000,1.763000);
positions[76] = Vec3(1.548000,1.522000,0.041000);
positions[77] = Vec3(1.485000,0.561000,0.369000);
positions[78] = Vec3(0.350000,1.661000,0.928000);
positions[79] = Vec3(1.653000,1.223000,0.578000);
positions[80] = Vec3(0.648000,1.349000,0.253000);
positions[81] = Vec3(0.340000,1.820000,0.483000);
positions[82] = Vec3(2.926000,0.119000,1.421000);
positions[83] = Vec3(1.512000,1.084000,0.156000);
positions[84] = Vec3(1.600000,2.115000,1.792000);
positions[85] = Vec3(1.089000,0.934000,1.584000);
positions[86] = Vec3(1.276000,1.104000,1.230000);
positions[87] = Vec3(0.485000,0.305000,0.428000);
positions[88] = Vec3(1.317000,1.261000,1.795000);
positions[89] = Vec3(0.039000,1.413000,1.085000);
positions[90] = Vec3(0.453000,0.701000,0.605000);
positions[91] = Vec3(1.283000,1.937000,0.752000);
positions[92] = Vec3(0.212000,1.416000,1.447000);
positions[93] = Vec3(0.203000,0.358000,0.723000);
positions[94] = Vec3(0.556000,0.445000,1.364000);
positions[95] = Vec3(1.436000,0.861000,0.911000);
positions[96] = Vec3(0.358000,0.966000,0.176000);
positions[97] = Vec3(1.478000,2.715000,0.427000);
positions[98] = Vec3(1.581000,0.575000,0.809000);
positions[99] = Vec3(1.007000,2.153000,2.887000);
positions[100] = Vec3(2.343000,0.663000,2.513000);
positions[101] = Vec3(2.105000,0.649000,1.635000);
positions[102] = Vec3(0.875000,0.743000,2.459000);
positions[103] = Vec3(0.229000,1.315000,1.879000);
positions[104] = Vec3(0.285000,0.935000,1.700000);
positions[105] = Vec3(2.269000,1.284000,2.234000);
positions[106] = Vec3(1.406000,1.149000,2.767000);
positions[107] = Vec3(1.076000,0.220000,1.849000);
positions[108] = Vec3(2.001000,1.532000,2.881000);
positions[109] = Vec3(2.893000,0.485000,1.860000);
positions[110] = Vec3(1.621000,1.786000,2.624000);
positions[111] = Vec3(0.500000,0.616000,1.818000);
positions[112] = Vec3(0.938000,2.978000,2.104000);
positions[113] = Vec3(0.550000,2.081000,0.454000);
positions[114] = Vec3(1.121000,0.685000,2.196000);
positions[115] = Vec3(1.088000,1.385000,2.184000);
positions[116] = Vec3(1.122000,2.705000,2.080000);
positions[117] = Vec3(0.918000,1.767000,2.861000);
positions[118] = Vec3(2.748000,0.232000,2.126000);
positions[119] = Vec3(1.238000,2.766000,0.109000);
positions[120] = Vec3(1.380000,0.785000,1.961000);
positions[121] = Vec3(1.236000,1.757000,0.150000);
positions[122] = Vec3(1.339000,2.187000,2.592000);
positions[123] = Vec3(1.414000,0.342000,2.714000);
positions[124] = Vec3(1.310000,0.770000,2.589000);
positions[125] = Vec3(1.686000,0.765000,2.321000);
positions[126] = Vec3(1.659000,1.367000,2.780000);
positions[127] = Vec3(0.141000,0.095000,1.903000);
positions[128] = Vec3(2.084000,1.002000,2.520000);
positions[129] = Vec3(2.819000,1.286000,2.626000);
positions[130] = Vec3(1.257000,1.044000,2.401000);
positions[131] = Vec3(1.064000,0.546000,2.839000);
positions[132] = Vec3(0.078000,1.246000,0.010000);
positions[133] = Vec3(1.506000,0.420000,2.223000);
positions[134] = Vec3(1.778000,0.699000,1.920000);
positions[135] = Vec3(1.315000,1.721000,2.733000);
positions[136] = Vec3(0.114000,0.281000,0.279000);
positions[137] = Vec3(1.082000,1.421000,2.596000);
positions[138] = Vec3(3.001000,0.592000,2.276000);
positions[139] = Vec3(1.384000,2.227000,2.992000);
positions[140] = Vec3(1.353000,0.044000,1.985000);
positions[141] = Vec3(1.367000,1.832000,2.383000);
positions[142] = Vec3(0.853000,1.119000,2.230000);
positions[143] = Vec3(1.675000,1.482000,2.295000);
positions[144] = Vec3(1.334000,1.890000,1.904000);
positions[145] = Vec3(1.630000,0.140000,2.939000);
positions[146] = Vec3(0.195000,1.290000,2.300000);
positions[147] = Vec3(2.178000,1.173000,3.001000);
positions[148] = Vec3(0.637000,0.655000,2.126000);
positions[149] = Vec3(0.993000,1.796000,2.529000);
positions[150] = Vec3(0.910000,0.701000,1.845000);
positions[151] = Vec3(0.191000,2.128000,0.355000);
positions[152] = Vec3(0.692000,1.163000,2.578000);
positions[153] = Vec3(1.154000,1.052000,1.974000);
positions[154] = Vec3(1.682000,0.335000,2.509000);
positions[155] = Vec3(0.569000,1.032000,1.895000);
positions[156] = Vec3(1.800000,2.796000,1.295000);
positions[157] = Vec3(2.517000,2.347000,2.878000);
positions[158] = Vec3(0.639000,2.470000,1.678000);
positions[159] = Vec3(0.634000,2.006000,1.829000);
positions[160] = Vec3(0.892000,0.215000,0.566000);
positions[161] = Vec3(1.800000,2.143000,1.491000);
positions[162] = Vec3(1.898000,0.226000,2.765000);
positions[163] = Vec3(0.791000,1.738000,0.260000);
positions[164] = Vec3(0.437000,1.740000,2.048000);
positions[165] = Vec3(1.687000,2.438000,1.166000);
positions[166] = Vec3(1.337000,2.304000,1.643000);
positions[167] = Vec3(1.270000,2.397000,1.033000);
positions[168] = Vec3(0.702000,2.429000,0.591000);
positions[169] = Vec3(0.842000,1.976000,0.724000);
positions[170] = Vec3(1.965000,0.095000,1.206000);
positions[171] = Vec3(0.355000,2.710000,0.618000);
positions[172] = Vec3(0.745000,1.434000,2.781000);
positions[173] = Vec3(0.707000,2.171000,1.502000);
positions[174] = Vec3(1.294000,2.696000,0.847000);
positions[175] = Vec3(1.143000,2.075000,0.276000);
positions[176] = Vec3(1.111000,2.474000,0.312000);
positions[177] = Vec3(1.144000,2.316000,0.633000);
positions[178] = Vec3(1.335000,0.292000,0.515000);
positions[179] = Vec3(1.926000,2.813000,2.703000);
positions[180] = Vec3(0.559000,2.314000,2.904000);
positions[181] = Vec3(1.308000,1.605000,1.534000);
positions[182] = Vec3(0.773000,2.913000,1.217000);
positions[183] = Vec3(1.612000,0.082000,1.027000);
positions[184] = Vec3(1.510000,0.287000,1.787000);
positions[185] = Vec3(0.716000,1.424000,1.843000);
positions[186] = Vec3(1.267000,0.685000,1.160000);
positions[187] = Vec3(0.306000,1.653000,1.717000);
positions[188] = Vec3(0.349000,0.020000,1.275000);
positions[189] = Vec3(0.166000,1.979000,0.804000);
positions[190] = Vec3(1.523000,2.992000,0.711000);
positions[191] = Vec3(1.998000,2.146000,0.088000);
positions[192] = Vec3(0.047000,2.513000,0.642000);
positions[193] = Vec3(0.501000,1.793000,1.438000);
positions[194] = Vec3(1.099000,2.010000,1.626000);
positions[195] = Vec3(2.580000,2.854000,1.328000);
positions[196] = Vec3(1.080000,2.779000,1.190000);
positions[197] = Vec3(0.901000,2.561000,0.948000);
positions[198] = Vec3(0.920000,2.990000,0.844000);
positions[199] = Vec3(0.819000,2.924000,1.711000);
positions[200] = Vec3(0.434000,1.516000,0.063000);
positions[201] = Vec3(1.470000,0.058000,0.231000);
positions[202] = Vec3(0.530000,3.005000,1.550000);
positions[203] = Vec3(0.447000,2.330000,1.277000);
positions[204] = Vec3(1.632000,2.683000,1.593000);
positions[205] = Vec3(0.885000,1.835000,2.072000);
positions[206] = Vec3(0.868000,2.601000,1.425000);
positions[207] = Vec3(0.720000,2.242000,0.907000);
positions[208] = Vec3(1.194000,0.144000,1.065000);
positions[209] = Vec3(0.448000,2.485000,0.959000);
positions[210] = Vec3(1.377000,2.694000,1.352000);
positions[211] = Vec3(1.305000,2.928000,2.713000);
positions[212] = Vec3(1.784000,2.456000,1.981000);
positions[213] = Vec3(0.354000,2.136000,1.563000);
positions[214] = Vec3(0.489000,2.000000,1.108000);
positions[215] = Vec3(1.884000,2.221000,0.461000);
positions[216] = Vec3(1.860000,2.540000,0.306000);
positions[217] = Vec3(1.753000,2.335000,2.768000);
positions[218] = Vec3(1.536000,2.441000,2.344000);
positions[219] = Vec3(0.531000,0.025000,2.235000);
positions[220] = Vec3(0.809000,0.011000,2.834000);
positions[221] = Vec3(0.289000,2.614000,2.879000);
positions[222] = Vec3(0.613000,1.891000,2.337000);
positions[223] = Vec3(0.507000,0.037000,2.694000);
positions[224] = Vec3(0.882000,2.185000,2.583000);
positions[225] = Vec3(0.503000,2.051000,2.615000);
positions[226] = Vec3(1.907000,1.956000,2.831000);
positions[227] = Vec3(2.833000,2.769000,1.644000);
positions[228] = Vec3(1.141000,0.113000,2.945000);
positions[229] = Vec3(0.600000,1.338000,2.200000);
positions[230] = Vec3(0.904000,2.360000,1.952000);
positions[231] = Vec3(0.738000,1.568000,2.437000);
positions[232] = Vec3(1.136000,2.535000,2.805000);
positions[233] = Vec3(1.430000,2.767000,2.321000);
positions[234] = Vec3(1.078000,2.470000,2.385000);
positions[235] = Vec3(0.296000,2.376000,2.560000);
positions[236] = Vec3(0.719000,0.300000,0.075000);
positions[237] = Vec3(0.518000,1.911000,0.080000);
positions[238] = Vec3(0.381000,1.570000,2.450000);
positions[239] = Vec3(0.716000,2.581000,2.697000);
positions[240] = Vec3(1.473000,2.617000,1.936000);
positions[241] = Vec3(0.421000,2.449000,0.229000);
positions[242] = Vec3(0.425000,2.817000,1.910000);
positions[243] = Vec3(1.312000,2.294000,2.057000);
positions[244] = Vec3(1.239000,0.007000,1.539000);
positions[245] = Vec3(0.822000,0.379000,2.086000);
positions[246] = Vec3(0.560000,2.562000,2.227000);
positions[247] = Vec3(0.863000,2.417000,0.050000);
positions[248] = Vec3(1.263000,0.151000,2.332000);
positions[249] = Vec3(0.237000,0.208000,2.336000);
positions[250] = Vec3(0.437000,2.370000,1.910000);
positions[251] = Vec3(1.119000,2.058000,2.207000);
positions[252] = Vec3(1.960000,1.749000,0.118000);
positions[253] = Vec3(2.415000,0.870000,2.757000);
positions[254] = Vec3(1.781000,0.342000,0.366000);
positions[255] = Vec3(2.172000,1.279000,1.421000);
positions[256] = Vec3(1.986000,0.715000,1.301000);
positions[257] = Vec3(1.657000,1.804000,0.810000);
positions[258] = Vec3(2.405000,1.202000,0.416000);
positions[259] = Vec3(1.932000,1.457000,0.786000);
positions[260] = Vec3(1.901000,1.271000,1.207000);
positions[261] = Vec3(1.864000,0.301000,0.810000);
positions[262] = Vec3(1.658000,0.673000,1.558000);
positions[263] = Vec3(2.637000,2.247000,0.396000);
positions[264] = Vec3(1.353000,0.369000,1.438000);
positions[265] = Vec3(0.530000,2.688000,1.346000);
positions[266] = Vec3(0.237000,0.485000,1.047000);
positions[267] = Vec3(2.806000,0.601000,0.822000);
positions[268] = Vec3(1.617000,2.018000,2.136000);
positions[269] = Vec3(2.000000,2.898000,0.022000);
positions[270] = Vec3(2.049000,1.883000,1.001000);
positions[271] = Vec3(2.477000,0.355000,1.786000);
positions[272] = Vec3(1.646000,0.983000,1.266000);
positions[273] = Vec3(1.683000,2.097000,1.114000);
positions[274] = Vec3(2.161000,0.921000,1.065000);
positions[275] = Vec3(2.099000,0.463000,1.942000);
positions[276] = Vec3(2.561000,1.638000,0.572000);
positions[277] = Vec3(2.205000,0.395000,1.005000);
positions[278] = Vec3(2.836000,0.203000,0.698000);
positions[279] = Vec3(2.662000,0.909000,0.966000);
positions[280] = Vec3(0.334000,0.350000,2.767000);
positions[281] = Vec3(2.241000,2.934000,1.248000);
positions[282] = Vec3(2.599000,2.953000,0.921000);
positions[283] = Vec3(2.219000,0.262000,0.058000);
positions[284] = Vec3(0.274000,0.656000,1.456000);
positions[285] = Vec3(1.814000,1.008000,0.882000);
positions[286] = Vec3(2.793000,1.395000,0.316000);
positions[287] = Vec3(0.773000,1.753000,1.639000);
positions[288] = Vec3(2.321000,0.994000,1.591000);
positions[289] = Vec3(2.243000,2.255000,1.690000);
positions[290] = Vec3(0.178000,1.342000,0.327000);
positions[291] = Vec3(1.623000,1.756000,1.426000);
positions[292] = Vec3(2.252000,0.109000,0.375000);
positions[293] = Vec3(3.003000,1.895000,1.895000);
positions[294] = Vec3(0.407000,0.831000,2.756000);
positions[295] = Vec3(2.193000,0.956000,0.632000);
positions[296] = Vec3(2.405000,0.641000,1.107000);
positions[297] = Vec3(2.361000,0.958000,0.162000);
positions[298] = Vec3(2.173000,1.544000,0.528000);
positions[299] = Vec3(1.565000,1.380000,1.428000);
positions[300] = Vec3(2.342000,0.538000,0.253000);
positions[301] = Vec3(1.910000,0.701000,0.954000);
positions[302] = Vec3(2.910000,0.288000,2.938000);
positions[303] = Vec3(0.257000,1.189000,0.958000);
positions[304] = Vec3(0.134000,1.775000,1.243000);
positions[305] = Vec3(2.476000,1.583000,1.956000);
positions[306] = Vec3(1.838000,1.791000,2.354000);
positions[307] = Vec3(1.906000,1.338000,1.696000);
positions[308] = Vec3(2.413000,2.869000,0.166000);
positions[309] = Vec3(3.006000,1.038000,1.322000);
positions[310] = Vec3(1.961000,0.962000,1.605000);
positions[311] = Vec3(0.082000,2.857000,0.020000);
positions[312] = Vec3(2.408000,1.499000,0.062000);
positions[313] = Vec3(2.349000,0.267000,1.415000);
positions[314] = Vec3(2.327000,1.717000,2.350000);
positions[315] = Vec3(2.928000,0.810000,1.582000);
positions[316] = Vec3(2.150000,0.336000,0.576000);
positions[317] = Vec3(2.664000,1.085000,2.962000);
positions[318] = Vec3(2.851000,0.670000,1.174000);
positions[319] = Vec3(1.954000,1.013000,1.975000);
positions[320] = Vec3(2.474000,1.542000,1.545000);
positions[321] = Vec3(2.826000,0.455000,1.490000);
positions[322] = Vec3(2.140000,2.826000,0.558000);
positions[323] = Vec3(2.151000,1.684000,1.780000);
positions[324] = Vec3(0.174000,0.673000,0.397000);
positions[325] = Vec3(0.066000,1.708000,0.160000);
positions[326] = Vec3(2.158000,0.303000,2.582000);
positions[327] = Vec3(2.602000,1.611000,2.632000);
positions[328] = Vec3(2.566000,1.138000,2.465000);
positions[329] = Vec3(2.026000,1.443000,2.477000);
positions[330] = Vec3(2.365000,0.309000,2.255000);
positions[331] = Vec3(1.636000,1.107000,2.058000);
positions[332] = Vec3(2.522000,2.584000,2.399000);
positions[333] = Vec3(2.537000,2.900000,2.158000);
positions[334] = Vec3(2.660000,0.537000,2.577000);
positions[335] = Vec3(2.679000,1.158000,1.724000);
positions[336] = Vec3(0.220000,1.894000,2.498000);
positions[337] = Vec3(2.266000,1.248000,1.837000);
positions[338] = Vec3(0.055000,1.656000,2.128000);
positions[339] = Vec3(2.899000,1.902000,2.823000);
positions[340] = Vec3(0.085000,2.994000,2.720000);
positions[341] = Vec3(0.013000,0.889000,2.468000);
positions[342] = Vec3(1.804000,0.372000,1.636000);
positions[343] = Vec3(0.201000,1.616000,2.824000);
positions[344] = Vec3(0.369000,1.273000,2.699000);
positions[345] = Vec3(2.996000,0.355000,2.596000);
positions[346] = Vec3(2.867000,1.314000,2.107000);
positions[347] = Vec3(2.611000,0.563000,2.140000);
positions[348] = Vec3(2.676000,2.954000,2.955000);
positions[349] = Vec3(0.256000,0.848000,2.062000);
positions[350] = Vec3(2.530000,0.028000,2.528000);
positions[351] = Vec3(0.537000,1.273000,1.596000);
positions[352] = Vec3(0.004000,1.004000,0.401000);
positions[353] = Vec3(1.676000,1.060000,2.463000);
positions[354] = Vec3(2.622000,1.473000,2.257000);
positions[355] = Vec3(2.917000,2.991000,2.316000);
positions[356] = Vec3(0.672000,1.123000,2.984000);
positions[357] = Vec3(2.229000,1.806000,2.673000);
positions[358] = Vec3(0.463000,0.951000,2.383000);
positions[359] = Vec3(2.126000,0.049000,2.037000);
positions[360] = Vec3(2.868000,0.876000,2.015000);
positions[361] = Vec3(2.720000,2.582000,0.079000);
positions[362] = Vec3(1.966000,0.693000,2.624000);
positions[363] = Vec3(1.971000,0.398000,2.318000);
positions[364] = Vec3(0.337000,0.630000,2.458000);
positions[365] = Vec3(2.562000,1.044000,2.040000);
positions[366] = Vec3(2.817000,1.485000,2.963000);
positions[367] = Vec3(2.514000,0.621000,2.992000);
positions[368] = Vec3(3.000000,1.551000,2.496000);
positions[369] = Vec3(0.698000,2.167000,2.180000);
positions[370] = Vec3(2.693000,0.849000,2.389000);
positions[371] = Vec3(2.092000,2.565000,2.986000);
positions[372] = Vec3(2.010000,3.001000,0.819000);
positions[373] = Vec3(2.392000,2.622000,1.636000);
positions[374] = Vec3(2.086000,2.325000,1.340000);
positions[375] = Vec3(2.578000,2.971000,0.502000);
positions[376] = Vec3(1.871000,2.789000,2.225000);
positions[377] = Vec3(2.230000,2.985000,1.594000);
positions[378] = Vec3(2.860000,2.788000,0.729000);
positions[379] = Vec3(2.051000,1.928000,1.472000);
positions[380] = Vec3(2.307000,2.219000,1.067000);
positions[381] = Vec3(2.369000,2.572000,1.289000);
positions[382] = Vec3(2.206000,1.924000,0.693000);
positions[383] = Vec3(1.984000,2.058000,2.005000);
positions[384] = Vec3(2.287000,1.854000,0.317000);
positions[385] = Vec3(2.525000,0.345000,0.686000);
positions[386] = Vec3(2.933000,1.920000,1.053000);
positions[387] = Vec3(0.324000,2.324000,0.601000);
positions[388] = Vec3(2.042000,1.576000,1.277000);
positions[389] = Vec3(0.031000,2.376000,0.949000);
positions[390] = Vec3(2.519000,2.250000,1.465000);
positions[391] = Vec3(0.221000,2.722000,1.652000);
positions[392] = Vec3(2.409000,2.361000,2.051000);
positions[393] = Vec3(2.472000,1.917000,1.673000);
positions[394] = Vec3(0.999000,2.715000,0.562000);
positions[395] = Vec3(1.669000,0.017000,1.508000);
positions[396] = Vec3(1.924000,1.777000,0.542000);
positions[397] = Vec3(2.635000,2.634000,1.905000);
positions[398] = Vec3(2.042000,2.628000,1.025000);
positions[399] = Vec3(2.694000,1.974000,2.009000);
positions[400] = Vec3(2.988000,2.221000,1.333000);
positions[401] = Vec3(1.772000,0.196000,1.978000);
positions[402] = Vec3(2.893000,2.961000,0.283000);
positions[403] = Vec3(2.615000,0.261000,0.245000);
positions[404] = Vec3(2.797000,2.521000,1.412000);
positions[405] = Vec3(0.013000,2.497000,0.246000);
positions[406] = Vec3(1.875000,2.861000,1.801000);
positions[407] = Vec3(2.800000,2.617000,1.049000);
positions[408] = Vec3(2.824000,1.858000,1.487000);
positions[409] = Vec3(2.434000,1.868000,1.275000);
positions[410] = Vec3(2.814000,0.526000,0.384000);
positions[411] = Vec3(2.844000,2.545000,2.246000);
positions[412] = Vec3(1.896000,2.587000,0.719000);
positions[413] = Vec3(0.350000,0.055000,0.076000);
positions[414] = Vec3(2.686000,1.784000,0.222000);
positions[415] = Vec3(2.724000,1.604000,0.989000);
positions[416] = Vec3(0.807000,1.761000,1.122000);
positions[417] = Vec3(2.120000,2.382000,2.226000);
positions[418] = Vec3(2.058000,1.587000,2.067000);
positions[419] = Vec3(2.904000,2.571000,2.686000);
positions[420] = Vec3(2.228000,2.910000,2.410000);
positions[421] = Vec3(2.797000,2.142000,0.114000);
positions[422] = Vec3(2.905000,1.875000,0.480000);
positions[423] = Vec3(1.881000,2.565000,2.469000);
positions[424] = Vec3(2.404000,1.929000,2.999000);
positions[425] = Vec3(2.389000,2.814000,2.782000);
positions[426] = Vec3(2.520000,0.301000,2.815000);
positions[427] = Vec3(2.726000,1.907000,2.339000);
positions[428] = Vec3(2.880000,2.273000,2.500000);
positions[429] = Vec3(2.574000,2.045000,2.716000);
positions[430] = Vec3(2.988000,2.288000,2.001000);
positions[431] = Vec3(0.011000,2.341000,2.904000);
positions[432] = Vec3(0.215000,2.265000,2.257000);
positions[433] = Vec3(2.268000,2.311000,0.234000);
positions[434] = Vec3(2.462000,2.621000,0.550000);
positions[435] = Vec3(1.530000,2.540000,2.728000);
positions[436] = Vec3(2.162000,2.306000,2.687000);
positions[437] = Vec3(2.748000,2.301000,1.734000);
positions[438] = Vec3(2.334000,1.976000,2.041000);
positions[439] = Vec3(1.981000,2.076000,2.443000);
positions[440] = Vec3(2.301000,1.367000,2.665000);
positions[441] = Vec3(2.399000,2.164000,2.403000);
positions[442] = Vec3(0.244000,2.713000,2.257000);
positions[443] = Vec3(0.683000,0.488000,2.781000);
positions[444] = Vec3(2.194000,2.711000,1.993000);
positions[445] = Vec3(2.947000,2.848000,2.001000);
positions[446] = Vec3(0.223000,1.981000,2.913000);
positions[447] = Vec3(0.010000,1.226000,0.917000);
positions[448] = Vec3(1.911000,0.426000,0.582000);
positions[449] = Vec3(2.204000,0.015000,0.136000);
positions[450] = Vec3(0.927000,0.138000,1.645000);
positions[451] = Vec3(0.155000,0.885000,1.479000);
positions[452] = Vec3(1.550000,1.933000,1.261000);
positions[453] = Vec3(1.304000,0.407000,0.287000);
positions[454] = Vec3(0.270000,1.384000,2.910000);
positions[455] = Vec3(0.516000,1.817000,1.695000);
positions[456] = Vec3(1.458000,2.879000,1.523000);
positions[457] = Vec3(1.702000,1.670000,0.593000);
positions[458] = Vec3(1.974000,1.380000,0.534000);
positions[459] = Vec3(2.835000,1.185000,0.479000);
positions[460] = Vec3(0.548000,2.979000,1.126000);
positions[461] = Vec3(1.202000,2.174000,1.466000);
positions[462] = Vec3(1.237000,1.701000,0.653000);
positions[463] = Vec3(2.939000,0.761000,0.349000);
positions[464] = Vec3(1.667000,2.119000,0.377000);
positions[465] = Vec3(1.196000,0.552000,1.372000);
positions[466] = Vec3(1.416000,0.901000,1.178000);
positions[467] = Vec3(0.519000,1.577000,2.227000);
positions[468] = Vec3(1.214000,1.281000,1.063000);
positions[469] = Vec3(0.822000,0.433000,1.375000);
positions[470] = Vec3(0.095000,2.760000,0.604000);
positions[471] = Vec3(1.325000,2.144000,1.848000);
positions[472] = Vec3(0.681000,0.896000,1.285000);
positions[473] = Vec3(0.406000,2.936000,0.717000);
positions[474] = Vec3(0.565000,1.852000,0.349000);
positions[475] = Vec3(0.597000,1.651000,1.020000);
positions[476] = Vec3(1.236000,0.170000,1.335000);
positions[477] = Vec3(0.586000,0.441000,1.980000);
positions[478] = Vec3(1.443000,1.208000,1.575000);
positions[479] = Vec3(0.247000,0.243000,0.502000);
positions[480] = Vec3(1.386000,1.564000,0.236000);
positions[481] = Vec3(0.871000,1.063000,0.930000);
positions[482] = Vec3(0.136000,0.992000,0.621000);
positions[483] = Vec3(0.889000,0.986000,0.010000);
positions[484] = Vec3(1.107000,2.731000,1.452000);
positions[485] = Vec3(0.942000,2.471000,0.517000);
positions[486] = Vec3(0.989000,0.652000,0.747000);
positions[487] = Vec3(0.899000,1.235000,2.707000);
positions[488] = Vec3(1.105000,0.684000,0.068000);
positions[489] = Vec3(1.660000,1.235000,2.276000);
positions[490] = Vec3(1.593000,1.883000,1.915000);
positions[491] = Vec3(1.528000,1.613000,0.920000);
positions[492] = Vec3(0.459000,1.046000,1.011000);
positions[493] = Vec3(0.213000,0.612000,0.644000);
positions[494] = Vec3(0.078000,1.392000,1.676000);
positions[495] = Vec3(0.605000,0.491000,0.574000);
positions[496] = Vec3(0.990000,1.586000,1.076000);
positions[497] = Vec3(0.297000,1.434000,1.028000);
positions[498] = Vec3(1.101000,1.471000,1.443000);
positions[499] = Vec3(0.072000,0.139000,1.653000);
positions[500] = Vec3(0.633000,0.884000,0.645000);
positions[501] = Vec3(0.352000,2.841000,1.463000);
positions[502] = Vec3(0.418000,0.774000,0.350000);
positions[503] = Vec3(2.641000,0.198000,0.869000);
positions[504] = Vec3(0.608000,1.341000,0.695000);
positions[505] = Vec3(1.778000,1.151000,1.830000);
positions[506] = Vec3(1.669000,0.342000,2.768000);
positions[507] = Vec3(1.256000,0.994000,0.798000);
positions[508] = Vec3(1.068000,0.375000,1.036000);
positions[509] = Vec3(0.910000,0.758000,1.589000);
positions[510] = Vec3(0.243000,2.452000,0.805000);
positions[511] = Vec3(1.018000,0.764000,1.122000);
positions[512] = Vec3(2.464000,1.089000,1.404000);
positions[513] = Vec3(0.670000,0.564000,0.034000);
positions[514] = Vec3(0.030000,1.296000,1.310000);
positions[515] = Vec3(1.210000,1.785000,1.691000);
positions[516] = Vec3(0.022000,0.620000,0.974000);
positions[517] = Vec3(1.499000,1.277000,2.986000);
positions[518] = Vec3(1.227000,1.896000,1.006000);
positions[519] = Vec3(0.528000,1.022000,1.635000);
positions[520] = Vec3(1.887000,2.670000,0.089000);
positions[521] = Vec3(1.661000,0.825000,0.793000);
positions[522] = Vec3(0.831000,1.494000,0.374000);
positions[523] = Vec3(1.356000,0.613000,0.930000);
positions[524] = Vec3(0.667000,0.600000,0.968000);
positions[525] = Vec3(1.154000,1.702000,2.925000);
positions[526] = Vec3(1.420000,1.581000,1.289000);
positions[527] = Vec3(1.383000,0.041000,0.932000);
positions[528] = Vec3(1.727000,0.140000,1.725000);
positions[529] = Vec3(0.711000,1.215000,2.004000);
positions[530] = Vec3(1.061000,1.067000,1.366000);
positions[531] = Vec3(0.377000,0.597000,1.224000);
positions[532] = Vec3(0.274000,0.719000,1.842000);
positions[533] = Vec3(0.840000,1.658000,1.874000);
positions[534] = Vec3(0.877000,0.290000,0.311000);
positions[535] = Vec3(2.130000,1.153000,1.196000);
positions[536] = Vec3(1.028000,1.379000,0.747000);
positions[537] = Vec3(1.107000,2.450000,2.079000);
positions[538] = Vec3(1.419000,1.333000,0.585000);
positions[539] = Vec3(0.430000,1.305000,1.369000);
positions[540] = Vec3(0.775000,1.363000,1.596000);
positions[541] = Vec3(1.522000,2.009000,0.736000);
positions[542] = Vec3(0.857000,1.722000,0.696000);
positions[543] = Vec3(0.722000,2.831000,1.478000);
positions[544] = Vec3(0.411000,1.673000,0.681000);
positions[545] = Vec3(1.511000,0.456000,0.597000);
positions[546] = Vec3(2.684000,0.820000,2.996000);
positions[547] = Vec3(1.593000,1.713000,2.369000);
positions[548] = Vec3(1.113000,0.803000,1.958000);
positions[549] = Vec3(1.267000,1.095000,0.254000);
positions[550] = Vec3(2.120000,0.540000,2.477000);
positions[551] = Vec3(0.566000,1.409000,2.588000);
positions[552] = Vec3(0.261000,0.872000,2.546000);
positions[553] = Vec3(1.878000,1.446000,2.680000);
positions[554] = Vec3(0.878000,1.606000,2.658000);
positions[555] = Vec3(1.564000,0.749000,1.786000);
positions[556] = Vec3(1.412000,1.942000,2.625000);
positions[557] = Vec3(1.660000,1.114000,2.710000);
positions[558] = Vec3(1.118000,0.813000,2.424000);
positions[559] = Vec3(1.482000,0.893000,2.434000);
positions[560] = Vec3(1.093000,1.129000,1.740000);
positions[561] = Vec3(2.163000,0.849000,2.709000);
positions[562] = Vec3(1.201000,1.429000,1.957000);
positions[563] = Vec3(0.235000,2.258000,2.002000);
positions[564] = Vec3(0.413000,1.444000,0.314000);
positions[565] = Vec3(0.164000,0.450000,2.408000);
positions[566] = Vec3(1.551000,0.851000,0.033000);
positions[567] = Vec3(0.659000,0.228000,2.807000);
positions[568] = Vec3(0.741000,0.131000,2.124000);
positions[569] = Vec3(0.455000,0.567000,2.682000);
positions[570] = Vec3(0.729000,0.971000,2.408000);
positions[571] = Vec3(1.487000,2.820000,0.162000);
positions[572] = Vec3(1.855000,0.700000,2.858000);
positions[573] = Vec3(0.305000,1.074000,1.926000);
positions[574] = Vec3(1.300000,0.153000,1.747000);
positions[575] = Vec3(1.272000,1.249000,2.568000);
positions[576] = Vec3(0.431000,0.743000,2.238000);
positions[577] = Vec3(0.493000,0.240000,0.184000);
positions[578] = Vec3(1.734000,0.506000,2.317000);
positions[579] = Vec3(0.874000,0.631000,2.692000);
positions[580] = Vec3(0.473000,2.790000,2.161000);
positions[581] = Vec3(1.310000,0.571000,2.759000);
positions[582] = Vec3(0.677000,0.798000,1.916000);
positions[583] = Vec3(0.944000,0.442000,1.858000);
positions[584] = Vec3(3.006000,2.098000,2.976000);
positions[585] = Vec3(0.864000,0.592000,2.231000);
positions[586] = Vec3(1.366000,0.611000,2.147000);
positions[587] = Vec3(2.871000,1.217000,2.880000);
positions[588] = Vec3(1.674000,2.664000,2.336000);
positions[589] = Vec3(1.757000,0.879000,2.101000);
positions[590] = Vec3(1.293000,2.939000,2.457000);
positions[591] = Vec3(1.108000,1.131000,2.206000);
positions[592] = Vec3(1.207000,1.658000,2.498000);
positions[593] = Vec3(0.850000,1.373000,2.312000);
positions[594] = Vec3(1.413000,1.060000,1.939000);
positions[595] = Vec3(1.138000,0.140000,2.102000);
positions[596] = Vec3(0.752000,1.307000,1.190000);
positions[597] = Vec3(1.254000,0.942000,2.790000);
positions[598] = Vec3(1.544000,1.614000,2.800000);
positions[599] = Vec3(2.128000,0.302000,2.833000);
positions[600] = Vec3(0.300000,1.744000,0.027000);
positions[601] = Vec3(1.878000,2.986000,2.060000);
positions[602] = Vec3(1.528000,0.233000,2.045000);
positions[603] = Vec3(1.146000,1.817000,2.067000);
positions[604] = Vec3(1.037000,2.746000,0.813000);
positions[605] = Vec3(0.524000,0.610000,1.566000);
positions[606] = Vec3(0.945000,2.964000,0.503000);
positions[607] = Vec3(1.788000,2.565000,0.965000);
positions[608] = Vec3(0.471000,2.510000,0.491000);
positions[609] = Vec3(0.512000,2.043000,1.371000);
positions[610] = Vec3(2.316000,2.423000,1.494000);
positions[611] = Vec3(1.575000,2.394000,2.953000);
positions[612] = Vec3(2.845000,2.869000,0.985000);
positions[613] = Vec3(1.016000,2.335000,1.003000);
positions[614] = Vec3(0.998000,2.830000,1.879000);
positions[615] = Vec3(0.624000,2.508000,0.075000);
positions[616] = Vec3(1.362000,2.808000,2.069000);
positions[617] = Vec3(1.747000,0.068000,0.810000);
positions[618] = Vec3(1.768000,2.355000,0.661000);
positions[619] = Vec3(1.535000,2.410000,2.085000);
positions[620] = Vec3(0.844000,2.004000,1.646000);
positions[621] = Vec3(1.124000,0.280000,0.649000);
positions[622] = Vec3(0.689000,2.170000,0.648000);
positions[623] = Vec3(0.849000,2.666000,1.175000);
positions[624] = Vec3(2.975000,1.963000,1.308000);
positions[625] = Vec3(1.074000,2.082000,0.714000);
positions[626] = Vec3(1.284000,2.651000,1.110000);
positions[627] = Vec3(1.669000,0.205000,0.180000);
positions[628] = Vec3(1.716000,0.047000,1.253000);
positions[629] = Vec3(0.501000,2.241000,1.043000);
positions[630] = Vec3(1.038000,1.833000,0.305000);
positions[631] = Vec3(0.646000,2.431000,1.424000);
positions[632] = Vec3(1.383000,2.059000,2.230000);
positions[633] = Vec3(0.370000,2.566000,1.192000);
positions[634] = Vec3(1.355000,2.006000,0.120000);
positions[635] = Vec3(2.113000,0.075000,0.589000);
positions[636] = Vec3(1.850000,0.448000,1.890000);
positions[637] = Vec3(1.215000,2.704000,0.405000);
positions[638] = Vec3(0.575000,2.997000,1.798000);
positions[639] = Vec3(0.967000,2.586000,2.603000);
positions[640] = Vec3(0.276000,1.669000,1.430000);
positions[641] = Vec3(1.483000,2.284000,1.128000);
positions[642] = Vec3(0.983000,3.003000,1.099000);
positions[643] = Vec3(0.539000,2.222000,1.720000);
positions[644] = Vec3(0.648000,2.826000,2.751000);
positions[645] = Vec3(0.803000,1.994000,0.993000);
positions[646] = Vec3(0.451000,0.216000,1.438000);
positions[647] = Vec3(1.604000,2.512000,0.334000);
positions[648] = Vec3(1.980000,2.022000,0.588000);
positions[649] = Vec3(1.843000,2.834000,1.544000);
positions[650] = Vec3(1.835000,3.005000,2.858000);
positions[651] = Vec3(0.679000,2.499000,0.838000);
positions[652] = Vec3(0.012000,2.637000,1.524000);
positions[653] = Vec3(0.314000,2.065000,0.602000);
positions[654] = Vec3(1.157000,0.004000,0.173000);
positions[655] = Vec3(0.736000,1.705000,1.382000);
positions[656] = Vec3(1.511000,2.736000,0.690000);
positions[657] = Vec3(1.330000,2.541000,1.735000);
positions[658] = Vec3(0.744000,0.170000,0.785000);
positions[659] = Vec3(2.593000,2.794000,0.703000);
positions[660] = Vec3(0.275000,1.872000,1.043000);
positions[661] = Vec3(1.624000,2.608000,1.341000);
positions[662] = Vec3(1.382000,0.122000,2.855000);
positions[663] = Vec3(1.326000,2.434000,0.783000);
positions[664] = Vec3(0.117000,0.116000,1.254000);
positions[665] = Vec3(1.045000,2.970000,2.748000);
positions[666] = Vec3(1.341000,2.692000,2.799000);
positions[667] = Vec3(1.797000,2.586000,2.709000);
positions[668] = Vec3(0.890000,2.484000,1.716000);
positions[669] = Vec3(2.373000,2.558000,1.889000);
positions[670] = Vec3(1.566000,2.323000,2.574000);
positions[671] = Vec3(1.257000,2.280000,0.399000);
positions[672] = Vec3(0.679000,2.130000,2.434000);
positions[673] = Vec3(2.016000,2.334000,2.462000);
positions[674] = Vec3(1.077000,2.213000,2.416000);
positions[675] = Vec3(0.581000,1.950000,2.081000);
positions[676] = Vec3(0.805000,2.315000,2.810000);
positions[677] = Vec3(0.844000,1.787000,2.322000);
positions[678] = Vec3(0.980000,2.205000,0.129000);
positions[679] = Vec3(2.468000,0.603000,2.740000);
positions[680] = Vec3(2.366000,2.403000,2.299000);
positions[681] = Vec3(0.337000,2.487000,2.329000);
positions[682] = Vec3(2.007000,2.793000,2.452000);
positions[683] = Vec3(1.072000,2.571000,0.063000);
positions[684] = Vec3(1.217000,2.283000,2.806000);
positions[685] = Vec3(0.459000,2.477000,2.728000);
positions[686] = Vec3(0.958000,1.975000,2.710000);
positions[687] = Vec3(0.914000,2.111000,2.052000);
positions[688] = Vec3(0.768000,2.958000,0.075000);
positions[689] = Vec3(0.474000,1.805000,2.533000);
positions[690] = Vec3(1.313000,2.552000,2.395000);
positions[691] = Vec3(1.853000,2.014000,2.229000);
positions[692] = Vec3(2.405000,2.230000,2.658000);
positions[693] = Vec3(0.727000,1.781000,0.016000);
positions[694] = Vec3(0.974000,2.791000,2.271000);
positions[695] = Vec3(0.438000,0.096000,2.457000);
positions[696] = Vec3(0.652000,2.392000,2.064000);
positions[697] = Vec3(1.972000,2.209000,2.834000);
positions[698] = Vec3(0.333000,0.141000,2.088000);
positions[699] = Vec3(1.813000,1.952000,0.063000);
positions[700] = Vec3(0.166000,2.838000,1.877000);
positions[701] = Vec3(1.772000,0.487000,0.951000);
positions[702] = Vec3(1.924000,1.404000,1.434000);
positions[703] = Vec3(2.734000,0.348000,1.712000);
positions[704] = Vec3(2.874000,0.729000,1.811000);
positions[705] = Vec3(1.841000,0.877000,1.137000);
positions[706] = Vec3(2.327000,1.491000,1.768000);
positions[707] = Vec3(1.916000,1.483000,1.057000);
positions[708] = Vec3(2.783000,0.850000,0.745000);
positions[709] = Vec3(1.829000,1.526000,0.085000);
positions[710] = Vec3(2.426000,1.082000,0.652000);
positions[711] = Vec3(1.645000,1.241000,1.217000);
positions[712] = Vec3(2.286000,0.725000,0.084000);
positions[713] = Vec3(2.755000,0.691000,1.421000);
positions[714] = Vec3(2.651000,0.591000,1.023000);
positions[715] = Vec3(2.040000,0.863000,0.442000);
positions[716] = Vec3(0.035000,0.109000,2.497000);
positions[717] = Vec3(0.127000,1.410000,0.572000);
positions[718] = Vec3(2.174000,0.357000,0.307000);
positions[719] = Vec3(2.705000,1.508000,0.758000);
positions[720] = Vec3(2.223000,1.407000,2.913000);
positions[721] = Vec3(2.528000,1.722000,1.088000);
positions[722] = Vec3(2.860000,0.345000,0.198000);
positions[723] = Vec3(2.580000,1.789000,1.479000);
positions[724] = Vec3(2.779000,0.295000,1.295000);
positions[725] = Vec3(0.097000,0.434000,2.826000);
positions[726] = Vec3(2.952000,1.654000,1.091000);
positions[727] = Vec3(0.119000,1.878000,0.343000);
positions[728] = Vec3(1.718000,1.173000,0.327000);
positions[729] = Vec3(2.833000,0.016000,0.527000);
positions[730] = Vec3(2.085000,1.779000,2.888000);
positions[731] = Vec3(2.754000,2.952000,1.485000);
positions[732] = Vec3(2.826000,0.935000,1.162000);
positions[733] = Vec3(1.564000,1.585000,1.615000);
positions[734] = Vec3(2.132000,0.645000,1.093000);
positions[735] = Vec3(2.294000,1.490000,1.350000);
positions[736] = Vec3(0.081000,0.490000,1.479000);
positions[737] = Vec3(2.118000,1.165000,1.642000);
positions[738] = Vec3(2.141000,0.121000,1.390000);
positions[739] = Vec3(2.385000,0.389000,1.196000);
positions[740] = Vec3(0.049000,0.166000,0.817000);
positions[741] = Vec3(1.993000,0.806000,1.814000);
positions[742] = Vec3(0.006000,1.450000,0.171000);
positions[743] = Vec3(2.297000,0.428000,0.764000);
positions[744] = Vec3(2.851000,0.469000,2.114000);
positions[745] = Vec3(1.814000,1.957000,0.945000);
positions[746] = Vec3(0.386000,0.327000,0.902000);
positions[747] = Vec3(2.452000,1.070000,1.807000);
positions[748] = Vec3(2.309000,1.537000,2.159000);
positions[749] = Vec3(2.712000,1.497000,2.007000);
positions[750] = Vec3(1.727000,0.924000,1.503000);
positions[751] = Vec3(0.861000,0.801000,0.344000);
positions[752] = Vec3(1.740000,1.245000,0.819000);
positions[753] = Vec3(0.117000,0.042000,0.197000);
positions[754] = Vec3(2.557000,0.996000,0.317000);
positions[755] = Vec3(2.228000,1.588000,2.548000);
positions[756] = Vec3(2.849000,1.557000,2.708000);
positions[757] = Vec3(0.152000,1.107000,0.219000);
positions[758] = Vec3(2.460000,1.318000,1.002000);
positions[759] = Vec3(2.405000,1.436000,0.528000);
positions[760] = Vec3(2.135000,1.179000,2.046000);
positions[761] = Vec3(1.726000,0.588000,0.286000);
positions[762] = Vec3(2.831000,1.053000,1.538000);
positions[763] = Vec3(1.932000,1.556000,1.833000);
positions[764] = Vec3(2.423000,0.900000,1.064000);
positions[765] = Vec3(3.001000,1.807000,0.709000);
positions[766] = Vec3(0.578000,1.095000,0.223000);
positions[767] = Vec3(2.215000,1.160000,0.252000);
positions[768] = Vec3(2.050000,0.921000,0.835000);
positions[769] = Vec3(2.080000,1.682000,0.738000);
positions[770] = Vec3(2.851000,1.753000,0.027000);
positions[771] = Vec3(0.203000,0.509000,0.202000);
positions[772] = Vec3(1.967000,1.018000,0.018000);
positions[773] = Vec3(1.869000,0.878000,2.472000);
positions[774] = Vec3(1.917000,0.228000,2.507000);
positions[775] = Vec3(0.316000,0.795000,2.991000);
positions[776] = Vec3(2.175000,1.229000,2.472000);
positions[777] = Vec3(2.405000,1.062000,2.931000);
positions[778] = Vec3(2.501000,0.511000,2.369000);
positions[779] = Vec3(2.641000,0.819000,2.141000);
positions[780] = Vec3(0.649000,1.384000,3.006000);
positions[781] = Vec3(1.012000,0.329000,2.963000);
positions[782] = Vec3(2.755000,0.350000,2.718000);
positions[783] = Vec3(2.315000,0.153000,2.454000);
positions[784] = Vec3(2.583000,1.696000,2.389000);
positions[785] = Vec3(0.439000,2.593000,1.776000);
positions[786] = Vec3(2.630000,1.390000,0.116000);
positions[787] = Vec3(2.854000,0.669000,2.478000);
positions[788] = Vec3(2.551000,1.342000,2.621000);
positions[789] = Vec3(2.533000,2.734000,2.987000);
positions[790] = Vec3(2.772000,2.446000,2.875000);
positions[791] = Vec3(2.817000,1.051000,2.498000);
positions[792] = Vec3(2.688000,1.404000,1.621000);
positions[793] = Vec3(0.083000,2.737000,2.775000);
positions[794] = Vec3(2.514000,0.322000,2.041000);
positions[795] = Vec3(2.470000,0.900000,2.504000);
positions[796] = Vec3(2.790000,0.444000,0.624000);
positions[797] = Vec3(0.040000,0.840000,2.202000);
positions[798] = Vec3(0.530000,1.067000,2.764000);
positions[799] = Vec3(0.191000,1.385000,2.541000);
positions[800] = Vec3(2.465000,0.363000,0.051000);
positions[801] = Vec3(1.850000,1.902000,2.592000);
positions[802] = Vec3(1.432000,0.306000,2.449000);
positions[803] = Vec3(2.259000,0.489000,1.753000);
positions[804] = Vec3(2.803000,1.118000,1.956000);
positions[805] = Vec3(2.426000,0.147000,1.636000);
positions[806] = Vec3(2.880000,1.846000,2.133000);
positions[807] = Vec3(2.862000,2.110000,1.867000);
positions[808] = Vec3(0.424000,1.184000,2.299000);
positions[809] = Vec3(2.518000,1.218000,2.228000);
positions[810] = Vec3(2.153000,0.834000,1.468000);
positions[811] = Vec3(0.105000,1.397000,2.088000);
positions[812] = Vec3(2.579000,0.601000,0.316000);
positions[813] = Vec3(2.594000,2.106000,2.968000);
positions[814] = Vec3(0.448000,1.435000,1.783000);
positions[815] = Vec3(2.125000,0.299000,2.132000);
positions[816] = Vec3(2.849000,1.402000,2.356000);
positions[817] = Vec3(2.956000,0.091000,2.078000);
positions[818] = Vec3(0.156000,1.696000,2.357000);
positions[819] = Vec3(1.566000,2.211000,1.557000);
positions[820] = Vec3(2.047000,0.194000,0.985000);
positions[821] = Vec3(1.947000,2.680000,0.488000);
positions[822] = Vec3(2.343000,2.796000,1.447000);
positions[823] = Vec3(2.006000,2.332000,0.265000);
positions[824] = Vec3(2.396000,1.834000,0.546000);
positions[825] = Vec3(2.538000,2.059000,2.207000);
positions[826] = Vec3(0.110000,2.360000,0.447000);
positions[827] = Vec3(2.198000,2.448000,1.136000);
positions[828] = Vec3(2.420000,2.121000,1.271000);
positions[829] = Vec3(0.422000,2.192000,0.260000);
positions[830] = Vec3(2.145000,2.767000,2.839000);
positions[831] = Vec3(2.434000,2.398000,0.421000);
positions[832] = Vec3(2.489000,2.175000,1.718000);
positions[833] = Vec3(2.870000,2.527000,0.814000);
positions[834] = Vec3(2.741000,2.016000,0.337000);
positions[835] = Vec3(1.997000,2.574000,2.107000);
positions[836] = Vec3(0.002000,2.128000,0.932000);
positions[837] = Vec3(2.787000,2.375000,0.234000);
positions[838] = Vec3(2.235000,1.852000,1.620000);
positions[839] = Vec3(2.782000,1.642000,0.422000);
positions[840] = Vec3(2.915000,1.760000,1.699000);
positions[841] = Vec3(2.047000,2.178000,1.549000);
positions[842] = Vec3(1.808000,1.878000,1.556000);
positions[843] = Vec3(2.224000,2.043000,0.913000);
positions[844] = Vec3(2.619000,2.611000,1.237000);
positions[845] = Vec3(2.916000,2.726000,0.168000);
positions[846] = Vec3(2.021000,2.833000,1.176000);
positions[847] = Vec3(2.967000,2.308000,2.258000);
positions[848] = Vec3(2.778000,2.270000,1.477000);
positions[849] = Vec3(2.121000,1.834000,2.002000);
positions[850] = Vec3(2.097000,2.752000,0.808000);
positions[851] = Vec3(1.897000,0.566000,1.501000);
positions[852] = Vec3(0.359000,2.802000,0.036000);
positions[853] = Vec3(2.966000,2.454000,1.186000);
positions[854] = Vec3(2.461000,2.964000,1.132000);
positions[855] = Vec3(2.093000,1.821000,1.243000);
positions[856] = Vec3(1.706000,2.659000,1.841000);
positions[857] = Vec3(2.074000,1.709000,0.342000);
positions[858] = Vec3(2.137000,2.894000,1.813000);
positions[859] = Vec3(0.223000,2.293000,1.417000);
positions[860] = Vec3(2.637000,0.007000,0.197000);
positions[861] = Vec3(1.416000,0.050000,0.483000);
positions[862] = Vec3(1.845000,2.250000,1.251000);
positions[863] = Vec3(2.906000,0.034000,2.896000);
positions[864] = Vec3(2.481000,0.204000,0.474000);
positions[865] = Vec3(2.234000,2.051000,0.158000);
positions[866] = Vec3(0.185000,2.453000,0.055000);
positions[867] = Vec3(2.509000,0.048000,2.786000);
positions[868] = Vec3(2.202000,2.206000,2.027000);
positions[869] = Vec3(0.061000,2.367000,2.656000);
positions[870] = Vec3(3.003000,2.755000,2.241000);
positions[871] = Vec3(0.297000,2.131000,2.463000);
positions[872] = Vec3(1.553000,0.429000,1.573000);
positions[873] = Vec3(2.506000,1.832000,1.911000);
positions[874] = Vec3(2.472000,1.814000,2.759000);
positions[875] = Vec3(1.922000,1.563000,2.278000);
positions[876] = Vec3(2.623000,2.666000,2.169000);
positions[877] = Vec3(0.120000,1.834000,2.723000);
positions[878] = Vec3(0.294000,0.103000,2.826000);
positions[879] = Vec3(2.364000,2.821000,0.417000);
positions[880] = Vec3(2.446000,1.734000,0.153000);
positions[881] = Vec3(2.777000,2.037000,2.565000);
positions[882] = Vec3(2.837000,2.477000,1.924000);
positions[883] = Vec3(2.221000,1.961000,2.443000);
positions[884] = Vec3(2.284000,2.895000,2.157000);
positions[885] = Vec3(2.728000,2.880000,1.861000);
positions[886] = Vec3(0.454000,2.080000,2.868000);
positions[887] = Vec3(2.430000,2.790000,2.524000);
positions[888] = Vec3(1.808000,2.213000,1.899000);
positions[889] = Vec3(2.666000,0.053000,2.309000);
positions[890] = Vec3(2.290000,2.408000,2.995000);
positions[891] = Vec3(2.646000,2.592000,1.625000);
positions[892] = Vec3(2.750000,2.508000,2.489000);
positions[893] = Vec3(0.211000,1.753000,1.939000);
platforms/opencl/tests/TestOpenCLCheckpoints.cpp
View file @ 2f553a66
......@@ -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-2013 Stanford University and the Authors. *
* Portions copyright (c) 2012-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,45 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests creating and loading checkpoints with the OpenCL platform.
*/
#include "OpenCLPlatform.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/AndersenThermostat.h"
#include "openmm/Context.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <sstream>
#include <vector>
using namespace OpenMM;
using namespace std;
static OpenCLPlatform platform;
const double TOL = 1e-5;
void compareStates(State& s1, State& s2) {
ASSERT_EQUAL_TOL(s1.getTime(), s2.getTime(), TOL);
int numParticles = s1.getPositions().size();
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(s1.getPositions()[i], s2.getPositions()[i], TOL);
ASSERT_EQUAL_VEC(s1.getVelocities()[i], s2.getVelocities()[i], TOL);
Vec3 a1, b1, c1, a2, b2, c2;
s1.getPeriodicBoxVectors(a1, b1, c1);
s2.getPeriodicBoxVectors(a2, b2, c2);
ASSERT_EQUAL_VEC(a1, a2, TOL);
ASSERT_EQUAL_VEC(b1, b2, TOL);
ASSERT_EQUAL_VEC(c1, c2, TOL);
for (map<string, double>::const_iterator iter = s1.getParameters().begin(); iter != s1.getParameters().end(); ++iter)
ASSERT_EQUAL(iter->second, (*s2.getParameters().find(iter->first)).second);
}
}
#include "OpenCLTests.h"
#include "TestCheckpoints.h"
void testCheckpoint() {
const int numParticles = 100;
......@@ -159,71 +122,6 @@ void testCheckpoint() {
compareStates(s6, s8);
}
void testSetState() {
const int numParticles = 10;
const double boxSize = 3.0;
const double temperature = 200.0;
System system;
system.addForce(new AndersenThermostat(0.0, 100.0));
NonbondedForce* nonbonded = new NonbondedForce();
system.addForce(nonbonded);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
nonbonded->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.2, 0.1);
positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
}
VerletIntegrator integrator(0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
context.setParameter(AndersenThermostat::Temperature(), temperature);
// Run for a little while.
integrator.step(100);
// Record the current state.
State s1 = context.getState(State::Positions | State::Velocities | State::Parameters);
// Continue the simulation for a few more steps and record a partial state.
integrator.step(10);
State s2 = context.getState(State::Positions);
// Restore the original state and see if everything gets restored correctly.
context.setPeriodicBoxVectors(Vec3(2*boxSize, 0, 0), Vec3(0, 2*boxSize, 0), Vec3(0, 0, 2*boxSize));
context.setParameter(AndersenThermostat::Temperature(), temperature+10);
context.setState(s1);
State s3 = context.getState(State::Positions | State::Velocities | State::Parameters);
compareStates(s1, s3);
// Set the partial state and see if the correct things were set.
context.setState(s2);
State s4 = context.getState(State::Positions | State::Velocities | State::Parameters);
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(s2.getPositions()[i], s4.getPositions()[i], TOL);
ASSERT_EQUAL_VEC(s1.getVelocities()[i], s4.getVelocities()[i], TOL);
}
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
testCheckpoint();
testSetState();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
void runPlatformTests() {
testCheckpoint();
}
platforms/opencl/tests/TestOpenCLEwald.cpp
View file @ 2f553a66
......@@ -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-2014 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,359 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the Ewald summation method OpenCL implementation of NonbondedForce.
*/
#include "OpenCLTests.h"
#include "TestEwald.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenCLPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/NonbondedForceImpl.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 testEwaldPME(bool includeExceptions) {
// Use amorphous NaCl system for the tests
const int numParticles = 894;
const double cutoff = 1.2;
const double boxSize = 3.00646;
double tol = 1e-5;
ReferencePlatform reference;
System system;
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(tol);
for (int i = 0; i < numParticles/2; i++)
system.addParticle(22.99);
for (int i = 0; i < numParticles/2; i++)
system.addParticle(35.45);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(1.0, 1.0,0.0);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(-1.0, 1.0,0.0);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
vector<Vec3> positions(numParticles);
#include "nacl_amorph.dat"
if (includeExceptions) {
// Add some exclusions.
for (int i = 0; i < numParticles-1; i++) {
Vec3 delta = positions[i]-positions[i+1];
if (sqrt(delta.dot(delta)) < 0.5*cutoff)
nonbonded->addException(i, i+1, i%2 == 0 ? 0.0 : 0.5, 1.0, 0.0);
}
}
// (1) Check whether the Reference and OpenCL platforms agree when using Ewald Method
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context clContext(system, integrator1, platform);
Context referenceContext(system, integrator2, reference);
clContext.setPositions(positions);
referenceContext.setPositions(positions);
State clState = clContext.getState(State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Forces | State::Energy);
tol = 1e-2;
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(referenceState.getForces()[i], clState.getForces()[i], tol);
}
tol = 1e-5;
ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), clState.getPotentialEnergy(), tol);
// (2) Check whether Ewald method in OpenCL is self-consistent
double norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = clState.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 5e-3;
double step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = clState.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
VerletIntegrator integrator3(0.01);
Context clContext2(system, integrator3, platform);
clContext2.setPositions(positions);
tol = 1e-2;
State clState2 = clContext2.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (clState2.getPotentialEnergy()-clState.getPotentialEnergy())/delta, tol)
// (3) Check whether the Reference and OpenCL platforms agree when using PME
nonbonded->setNonbondedMethod(NonbondedForce::PME);
clContext.reinitialize();
referenceContext.reinitialize();
clContext.setPositions(positions);
referenceContext.setPositions(positions);
clState = clContext.getState(State::Forces | State::Energy);
referenceState = referenceContext.getState(State::Forces | State::Energy);
tol = 1e-2;
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(referenceState.getForces()[i], clState.getForces()[i], tol);
}
tol = 1e-5;
ASSERT_EQUAL_TOL(referenceState.getPotentialEnergy(), clState.getPotentialEnergy(), tol);
// (4) Check whether PME method in OpenCL is self-consistent
norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = clState.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = clState.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
VerletIntegrator integrator4(0.01);
Context clContext3(system, integrator4, platform);
clContext3.setPositions(positions);
tol = 1e-2;
State clState3 = clContext3.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (clState3.getPotentialEnergy()-clState.getPotentialEnergy())/delta, tol)
}
void testEwald2Ions() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(-1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(TOL);
system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(3.048000,2.764000,3.156000);
positions[1] = Vec3(2.809000,2.888000,2.571000);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
ASSERT_EQUAL_VEC(Vec3( 123.711, -64.1877, 302.716), forces[1], 10*TOL);
ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 0.01/*10*TOL*/);
void runPlatformTests() {
}
void testTriclinic() {
// Create a triclinic box containing eight particles.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(2.5, 0, 0), Vec3(0.5, 3.0, 0), Vec3(0.7, 0.9, 3.5));
for (int i = 0; i < 8; i++)
system.addParticle(1.0);
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
force->setNonbondedMethod(NonbondedForce::PME);
force->setCutoffDistance(1.0);
force->setPMEParameters(3.45891, 32, 40, 48);
for (int i = 0; i < 4; i++)
force->addParticle(-1, 0.440104, 0.4184); // Cl parameters
for (int i = 0; i < 4; i++)
force->addParticle(1, 0.332840, 0.0115897); // Na parameters
vector<Vec3> positions(8);
positions[0] = Vec3(1.744, 2.788, 3.162);
positions[1] = Vec3(1.048, 0.762, 2.340);
positions[2] = Vec3(2.489, 1.570, 2.817);
positions[3] = Vec3(1.027, 1.893, 3.271);
positions[4] = Vec3(0.937, 0.825, 0.009);
positions[5] = Vec3(2.290, 1.887, 3.352);
positions[6] = Vec3(1.266, 1.111, 2.894);
positions[7] = Vec3(0.933, 1.862, 3.490);
// Compute the forces and energy.
VerletIntegrator integ(0.001);
Context context(system, integ, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
// Compare them to values computed by Gromacs.
double expectedEnergy = -963.370;
vector<Vec3> expectedForce(8);
expectedForce[0] = Vec3(4.25253e+01, -1.23503e+02, 1.22139e+02);
expectedForce[1] = Vec3(9.74752e+01, 1.68213e+02, 1.93169e+02);
expectedForce[2] = Vec3(-1.50348e+02, 1.29165e+02, 3.70435e+02);
expectedForce[3] = Vec3(9.18644e+02, -3.52571e+00, -1.34772e+03);
expectedForce[4] = Vec3(-1.61193e+02, 9.01528e+01, -7.12904e+01);
expectedForce[5] = Vec3(2.82630e+02, 2.78029e+01, -3.72864e+02);
expectedForce[6] = Vec3(-1.47454e+02, -2.14448e+02, -3.55789e+02);
expectedForce[7] = Vec3(-8.82195e+02, -7.39132e+01, 1.46202e+03);
for (int i = 0; i < 8; i++) {
ASSERT_EQUAL_VEC(expectedForce[i], state.getForces()[i], 1e-4);
}
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-4);
}
void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
// Create a cloud of random point charges.
const int numParticles = 51;
const double boxWidth = 5.0;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
}
force->setNonbondedMethod(method);
// For various values of the cutoff and error tolerance, see if the actual error is reasonable.
for (double cutoff = 1.0; cutoff < boxWidth/2; cutoff *= 1.2) {
force->setCutoffDistance(cutoff);
vector<Vec3> refForces;
double norm = 0.0;
for (double tol = 5e-5; tol < 1e-3; tol *= 2.0) {
force->setEwaldErrorTolerance(tol);
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces);
if (refForces.size() == 0) {
refForces = state.getForces();
for (int i = 0; i < numParticles; i++)
norm += refForces[i].dot(refForces[i]);
norm = sqrt(norm);
}
else {
double diff = 0.0;
for (int i = 0; i < numParticles; i++) {
Vec3 delta = refForces[i]-state.getForces()[i];
diff += delta.dot(delta);
}
diff = sqrt(diff)/norm;
ASSERT(diff < 2*tol);
}
if (method == NonbondedForce::PME) {
// See if the PME parameters were calculated correctly.
double expectedAlpha, actualAlpha;
int expectedSize[3], actualSize[3];
NonbondedForceImpl::calcPMEParameters(system, *force, expectedAlpha, expectedSize[0], expectedSize[1], expectedSize[2]);
force->getPMEParametersInContext(context, actualAlpha, actualSize[0], actualSize[1], actualSize[2]);
ASSERT_EQUAL_TOL(expectedAlpha, actualAlpha, 1e-5);
for (int i = 0; i < 3; i++) {
ASSERT(actualSize[i] >= expectedSize[i]);
ASSERT(actualSize[i] < expectedSize[i]+10);
}
}
}
}
}
void testPMEParameters() {
// Create a cloud of random point charges.
const int numParticles = 51;
const double boxWidth = 4.7;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
NonbondedForce* force = new NonbondedForce();
system.addForce(force);
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
}
force->setNonbondedMethod(NonbondedForce::PME);
// Compute the energy with an error tolerance of 1e-3.
force->setEwaldErrorTolerance(1e-3);
VerletIntegrator integrator1(0.01);
Context context1(system, integrator1, platform);
context1.setPositions(positions);
double energy1 = context1.getState(State::Energy).getPotentialEnergy();
// Try again with an error tolerance of 1e-4.
force->setEwaldErrorTolerance(1e-4);
VerletIntegrator integrator2(0.01);
Context context2(system, integrator2, platform);
context2.setPositions(positions);
double energy2 = context2.getState(State::Energy).getPotentialEnergy();
// Now explicitly set the parameters. These should match the values that were
// used for tolerance 1e-3.
force->setPMEParameters(2.49291157051793, 32, 32, 32);
VerletIntegrator integrator3(0.01);
Context context3(system, integrator3, platform);
context3.setPositions(positions);
double energy3 = context3.getState(State::Energy).getPotentialEnergy();
ASSERT_EQUAL_TOL(energy1, energy3, 1e-6);
ASSERT(fabs((energy1-energy2)/energy1) > 1e-5);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
testEwaldPME(false);
testEwaldPME(true);
// testEwald2Ions();
testTriclinic();
testErrorTolerance(NonbondedForce::Ewald);
testErrorTolerance(NonbondedForce::PME);
testPMEParameters();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
platforms/opencl/tests/TestOpenCLGBSAOBCForce.cpp
View file @ 2f553a66
......@@ -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-2014 Stanford University and the Authors. *
* Portions copyright (c) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -29,243 +29,8 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the OpenCL implementation of GBSAOBCForce.
*/
#include "OpenCLTests.h"
#include "TestGBSAOBCForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenCLPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/GBSAOBCForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include "openmm/NonbondedForce.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
static OpenCLPlatform platform;
const double TOL = 1e-5;
void testSingleParticle() {
System system;
system.addParticle(2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* gbsa = new GBSAOBCForce();
NonbondedForce* nonbonded = new NonbondedForce();
gbsa->addParticle( 0.5, 0.15, 1);
nonbonded->addParticle(0.5, 1, 0);
system.addForce(gbsa);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double bornRadius = 0.15-0.009; // dielectric offset
double eps0 = EPSILON0;
double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
double extendedRadius = 0.15+0.14; // probe radius
double nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
// Change the parameters and see if it is still correct.
gbsa->setParticleParameters(0, 0.4, 0.25, 1);
gbsa->updateParametersInContext(context);
state = context.getState(State::Energy);
bornRadius = 0.25-0.009; // dielectric offset
bornEnergy = (-0.4*0.4/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
extendedRadius = 0.25+0.14;
nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.25/bornRadius, 6.0);
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
}
void testGlobalSettings() {
System system;
system.addParticle(2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* forceField = new GBSAOBCForce();
forceField->addParticle(0.5, 0.15, 1);
const double soluteDielectric = 2.1;
const double solventDielectric = 35.0;
const double surfaceAreaEnergy = 0.75;
forceField->setSoluteDielectric(soluteDielectric);
forceField->setSolventDielectric(solventDielectric);
forceField->setSurfaceAreaEnergy(surfaceAreaEnergy);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double bornRadius = 0.15-0.009; // dielectric offset
double eps0 = EPSILON0;
double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/soluteDielectric-1.0/solventDielectric)/bornRadius;
double extendedRadius = 0.15+0.14; // probe radius
double nonpolarEnergy = 4*PI_M*surfaceAreaEnergy*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
}
void testCutoffAndPeriodic() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* gbsa = new GBSAOBCForce();
NonbondedForce* nonbonded = new NonbondedForce();
gbsa->addParticle(-1, 0.15, 1);
nonbonded->addParticle(-1, 1, 0);
gbsa->addParticle(1, 0.15, 1);
nonbonded->addParticle(1, 1, 0);
const double cutoffDistance = 3.0;
const double boxSize = 10.0;
nonbonded->setCutoffDistance(cutoffDistance);
gbsa->setCutoffDistance(cutoffDistance);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(gbsa);
system.addForce(nonbonded);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
// Calculate the forces for both cutoff and periodic with two different atom positions.
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
Context context(system, integrator, platform);
context.setPositions(positions);
State state1 = context.getState(State::Forces);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
context.reinitialize();
context.setPositions(positions);
State state2 = context.getState(State::Forces);
positions[1][0]+= boxSize;
nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffNonPeriodic);
context.reinitialize();
context.setPositions(positions);
State state3 = context.getState(State::Forces);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
gbsa->setNonbondedMethod(GBSAOBCForce::CutoffPeriodic);
context.reinitialize();
context.setPositions(positions);
State state4 = context.getState(State::Forces);
// All forces should be identical, exception state3 which should be zero.
ASSERT_EQUAL_VEC(state1.getForces()[0], state2.getForces()[0], 0.01);
ASSERT_EQUAL_VEC(state1.getForces()[1], state2.getForces()[1], 0.01);
ASSERT_EQUAL_VEC(state1.getForces()[0], state4.getForces()[0], 0.01);
ASSERT_EQUAL_VEC(state1.getForces()[1], state4.getForces()[1], 0.01);
ASSERT_EQUAL_VEC(state3.getForces()[0], Vec3(0, 0, 0), 0.01);
ASSERT_EQUAL_VEC(state3.getForces()[1], Vec3(0, 0, 0), 0.01);
void runPlatformTests() {
}
void testForce(int numParticles, NonbondedForce::NonbondedMethod method, GBSAOBCForce::NonbondedMethod method2) {
ReferencePlatform reference;
System system;
GBSAOBCForce* gbsa = new GBSAOBCForce();
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(1.0);
double charge = i%2 == 0 ? -1 : 1;
gbsa->addParticle(charge, 0.15, 1);
nonbonded->addParticle(charge, 1, 0);
}
nonbonded->setNonbondedMethod(method);
gbsa->setNonbondedMethod(method2);
nonbonded->setCutoffDistance(3.0);
gbsa->setCutoffDistance(3.0);
int grid = (int) floor(0.5+pow(numParticles, 1.0/3.0));
if (method == NonbondedForce::CutoffPeriodic) {
double boxSize = (grid+1)*1.1;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
}
system.addForce(gbsa);
system.addForce(nonbonded);
LangevinIntegrator integrator1(0, 0.1, 0.01);
LangevinIntegrator integrator2(0, 0.1, 0.01);
Context context(system, integrator1, platform);
Context refContext(system, integrator2, reference);
// Set random (but uniformly distributed) positions for all the particles.
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < grid; i++)
for (int j = 0; j < grid; j++)
for (int k = 0; k < grid; k++)
positions[i*grid*grid+j*grid+k] = Vec3(i*1.1, j*1.1, k*1.1);
for (int i = 0; i < numParticles; ++i)
positions[i] = positions[i] + Vec3(0.5*genrand_real2(sfmt), 0.5*genrand_real2(sfmt), 0.5*genrand_real2(sfmt));
context.setPositions(positions);
refContext.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
State refState = refContext.getState(State::Forces | State::Energy);
// Make sure the OpenCL and Reference platforms agree.
double norm = 0.0;
double diff = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = state.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
Vec3 delta = f-refState.getForces()[i];
diff += delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
}
norm = std::sqrt(norm);
diff = std::sqrt(diff);
ASSERT_EQUAL_TOL(0.0, diff, 0.001*norm);
ASSERT_EQUAL_TOL(state.getPotentialEnergy(), refState.getPotentialEnergy(), 1e-3);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
// (This doesn't work with cutoffs, since the energy changes discontinuously at the cutoff distance.)
if (method == NonbondedForce::NoCutoff)
{
const double delta = 0.3;
double step = 0.5*delta/norm;
vector<Vec3> positions2(numParticles), positions3(numParticles);
for (int i = 0; i < numParticles; ++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);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/delta, 1e-2)
}
}
int main() {
try {
testSingleParticle();
testGlobalSettings();
testCutoffAndPeriodic();
for (int i = 5; i < 11; i++) {
testForce(i*i*i, NonbondedForce::NoCutoff, GBSAOBCForce::NoCutoff);
testForce(i*i*i, NonbondedForce::CutoffNonPeriodic, GBSAOBCForce::CutoffNonPeriodic);
testForce(i*i*i, NonbondedForce::CutoffPeriodic, GBSAOBCForce::CutoffPeriodic);
}
}
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