Skip to content

GitLab

Commit 30fa6ecb authored by Peter Eastman's avatar Peter Eastman
Browse files

Improved test cases for variable time step integrators

parent 042b8c2c
Show whitespace changes
Inline Side-by-side
Showing with 362 additions and 14 deletions
platforms/cuda/tests/TestCudaVariableLangevinIntegrator.cpp
View file @ 30fa6ecb
......@@ -99,7 +99,7 @@ void testTemperature() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 2.0, 1e-4);
VariableLangevinIntegrator integrator(temp, 5.0, 1e-4);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
......@@ -122,11 +122,11 @@ void testTemperature() {
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
integrator.step(5);
}
ke /= 1000;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.1);
}
void testConstraints() {
......@@ -228,6 +228,65 @@ void testRandomSeed() {
}
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
int numParticles = 0;
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
++numParticles;
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableLangevinIntegrator integrator(temp, 6.0, 1e-5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(2.0);
// Make sure the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 400; ++i) {
double t = 2.0 + 0.01 * (i + 1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
}
ke /= 400;
double expected = 1.5 * numParticles * BOLTZ * temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.01);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
......@@ -236,6 +295,7 @@ int main(int argc, char* argv[]) {
testTemperature();
testConstraints();
testRandomSeed();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/cuda/tests/TestCudaVariableVerletIntegrator.cpp
View file @ 30fa6ecb
......@@ -211,6 +211,62 @@ void testConstrainedClusters() {
ASSERT(context.getState(State::Positions).getTime() > 0.1);
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableVerletIntegrator integrator(1e-5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(1.0);
// Simulate it and see whether energy remains constant.
State state0 = context.getState(State::Energy);
double initialEnergy = state0.getKineticEnergy() + state0.getPotentialEnergy();
for (int i = 0; i < 20; i++) {
double t = 1.0 + 0.05*(i+1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
double energy = state.getKineticEnergy() + state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
}
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
......@@ -218,6 +274,7 @@ int main(int argc, char* argv[]) {
testSingleBond();
testConstraints();
testConstrainedClusters();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/opencl/tests/TestOpenCLVariableLangevinIntegrator.cpp
View file @ 30fa6ecb
......@@ -99,7 +99,7 @@ void testTemperature() {
const int numParticles = 8;
const double temp = 100.0;
System system;
VariableLangevinIntegrator integrator(temp, 2.0, 1e-4);
VariableLangevinIntegrator integrator(temp, 5.0, 1e-4);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
......@@ -122,11 +122,11 @@ void testTemperature() {
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
integrator.step(5);
}
ke /= 1000;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.1);
}
void testConstraints() {
......@@ -228,6 +228,65 @@ void testRandomSeed() {
}
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
int numParticles = 0;
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
++numParticles;
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableLangevinIntegrator integrator(temp, 6.0, 1e-5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(2.0);
// Make sure the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 400; ++i) {
double t = 2.0 + 0.01 * (i + 1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
}
ke /= 400;
double expected = 1.5 * numParticles * BOLTZ * temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.01);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
......@@ -236,6 +295,7 @@ int main(int argc, char* argv[]) {
testTemperature();
testConstraints();
testRandomSeed();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......@@ -244,5 +304,3 @@ int main(int argc, char* argv[]) {
cout << "Done" << endl;
return 0;
}
platforms/opencl/tests/TestOpenCLVariableVerletIntegrator.cpp
View file @ 30fa6ecb
......@@ -210,6 +210,62 @@ void testConstrainedClusters() {
ASSERT(context.getState(State::Positions).getTime() > 0.1);
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableVerletIntegrator integrator(1e-5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(1.0);
// Simulate it and see whether energy remains constant.
State state0 = context.getState(State::Energy);
double initialEnergy = state0.getKineticEnergy() + state0.getPotentialEnergy();
for (int i = 0; i < 20; i++) {
double t = 1.0 + 0.05*(i+1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
double energy = state.getKineticEnergy() + state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
}
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
......@@ -217,6 +273,7 @@ int main(int argc, char* argv[]) {
testSingleBond();
testConstraints();
testConstrainedClusters();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......@@ -225,5 +282,3 @@ int main(int argc, char* argv[]) {
cout << "Done" << endl;
return 0;
}
platforms/reference/tests/TestReferenceVariableLangevinIntegrator.cpp
View file @ 30fa6ecb
......@@ -99,7 +99,7 @@ void testTemperature() {
const double temp = 100.0;
ReferencePlatform platform;
System system;
VariableLangevinIntegrator integrator(temp, 2.0, 1e-4);
VariableLangevinIntegrator integrator(temp, 5.0, 1e-4);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(2.0);
......@@ -122,11 +122,11 @@ void testTemperature() {
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
integrator.step(5);
}
ke /= 1000;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.1);
}
void testConstraints() {
......@@ -230,12 +230,73 @@ void testRandomSeed() {
}
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
ReferencePlatform platform;
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
int numParticles = 0;
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
++numParticles;
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableLangevinIntegrator integrator(temp, 6.0, 1e-5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(2.0);
// Make sure the temperature is correct.
double ke = 0.0;
for (int i = 0; i < 400; ++i) {
double t = 2.0 + 0.01 * (i + 1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
}
ke /= 400;
double expected = 1.5 * numParticles * BOLTZ * temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 0.01);
}
int main() {
try {
testSingleBond();
testTemperature();
testConstraints();
testRandomSeed();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......@@ -244,4 +305,3 @@ int main() {
cout << "Done" << endl;
return 0;
}
platforms/reference/tests/TestReferenceVariableVerletIntegrator.cpp
View file @ 30fa6ecb
......@@ -210,11 +210,69 @@ void testConstrainedClusters() {
ASSERT(context.getState(State::Positions).getTime() > 0.1);
}
void testArgonBox() {
const int gridSize = 8;
const double mass = 40.0; // Ar atomic mass
const double temp = 120.0; // K
const double epsilon = BOLTZ * temp; // L-J well depth for Ar
const double sigma = 0.34; // L-J size for Ar in nm
const double density = 0.8; // atoms / sigma^3
double cellSize = sigma / pow(density, 0.333);
double boxSize = gridSize * cellSize;
double cutoff = 2.0 * sigma;
// Create a box of argon atoms.
ReferencePlatform platform;
System system;
NonbondedForce* nonbonded = new NonbondedForce();
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
const Vec3 half(0.5, 0.5, 0.5);
for (int i = 0; i < gridSize; i++) {
for (int j = 0; j < gridSize; j++) {
for (int k = 0; k < gridSize; k++) {
system.addParticle(mass);
nonbonded->addParticle(0, sigma, epsilon);
positions.push_back((Vec3(i, j, k) + half + Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt))*0.1) * cellSize);
}
}
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
VariableVerletIntegrator integrator(1e-5);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temp);
// Equilibrate.
integrator.stepTo(1.0);
// Simulate it and see whether energy remains constant.
State state0 = context.getState(State::Energy);
double initialEnergy = state0.getKineticEnergy() + state0.getPotentialEnergy();
for (int i = 0; i < 20; i++) {
double t = 1.0 + 0.05*(i+1);
integrator.stepTo(t);
State state = context.getState(State::Energy);
double energy = state.getKineticEnergy() + state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
}
}
int main() {
try {
testSingleBond();
testConstraints();
testConstrainedClusters();
testArgonBox();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment