Skip to content

GitLab

Commit 49d60231 authored by Lee-Ping Wang's avatar Lee-Ping Wang
Browse files

Merge branch 'master' of https://github.com/leeping/openmm

parents 5587fac9 3a2f4643
Show whitespace changes
Inline Side-by-side
Showing with 509 additions and 391 deletions
openmmapi/include/openmm/MonteCarloAnisotropicBarostat.h
View file @ 49d60231
......@@ -9,8 +9,8 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2010 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Portions copyright (c) 2010-2013 Stanford University and the Authors. *
* Authors: Peter Eastman, Lee-Ping Wang *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
......@@ -43,8 +43,11 @@ namespace OpenMM {
* This class uses a Monte Carlo algorithm to adjust the size of the periodic box, simulating the
* effect of constant pressure.
*
* Compared to MonteCarloBarostat, this class scales the three axes of the simulation cell independently.
* The user supplies three doubles to specify the pressure along each axis.
* This class is similar to MonteCarloBarostat, but each Monte Carlo move is applied to only one axis
* of the periodic box (unlike MonteCarloBarostat, which scales the entire box isotropically). This
* means that the box may change shape as well as size over the course of the simulation. It also
* allows you to specify a different pressure for each axis of the box, or to keep the box size fixed
* along certain axes while still allowing it to change along others.
*
* This class assumes the simulation is also being run at constant temperature, and requires you
* to specify the system temperature (since it affects the acceptance probability for Monte Carlo
......@@ -84,39 +87,33 @@ public:
* @param defaultPressure The default pressure acting on each axis (in bar)
* @param temperature the temperature at which the system is being maintained (in Kelvin)
* @param frequency the frequency at which Monte Carlo pressure changes should be attempted (in time steps)
* @param scaleX on/off switch for whether to scale the X axis
* @param scaleY on/off switch for whether to scale the Y axis
* @param scaleZ on/off switch for whether to scale the Z axis
* @param scaleX whether to allow the X dimension of the periodic box to change size
* @param scaleY whether to allow the Y dimension of the periodic box to change size
* @param scaleZ whether to allow the Z dimension of the periodic box to change size
*/
MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double temperature, int frequency = 25, bool scaleX = 1, bool scaleY = 1, bool scaleZ = 1);
/**
* Get the default pressure (in bar).
*
* @return the default pressure acting on the system, measured in bar.
* @return the default pressure acting along each axis, measured in bar.
*/
Vec3 getDefaultPressure() const {
const Vec3& getDefaultPressure() const {
return defaultPressure;
}
/**
* Get the true/false flag for scaling the X-axis.
*
* @return the true/false flag for scaling the X-axis.
* Get whether to allow the X dimension of the periodic box to change size.
*/
bool getScaleX() const {
return scaleX;
}
/**
* Get the true/false flag for scaling the Y-axis.
*
* @return the true/false flag for scaling the Y-axis.
* Get whether to allow the Y dimension of the periodic box to change size.
*/
bool getScaleY() const {
return scaleY;
}
/**
* Get the true/false flag for scaling the Z-axis.
*
* @return the true/false flag for scaling the Z-axis.
* Get whether to allow the Z dimension of the periodic box to change size.
*/
bool getScaleZ() const {
return scaleZ;
......@@ -172,15 +169,7 @@ private:
double temperature;
bool scaleX, scaleY, scaleZ;
int frequency, randomNumberSeed;
double GetTemperature() const {
return temperature;
}
void SetTemperature(double temperature) {
this->temperature = temperature;
}
};
};
} // namespace OpenMM
......
openmmapi/include/openmm/MonteCarloBarostat.h
View file @ 49d60231
......@@ -123,15 +123,7 @@ protected:
private:
double defaultPressure, temperature;
int frequency, randomNumberSeed;
double GetTemperature() const {
return temperature;
}
void SetTemperature(double temperature) {
this->temperature = temperature;
}
};
};
} // namespace OpenMM
......
openmmapi/include/openmm/internal/AssertionUtilities.h
View file @ 49d60231
......@@ -9,7 +9,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) 2008-2013 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -54,6 +54,8 @@ void OPENMM_EXPORT throwException(const char* file, int line, const std::string&
#define ASSERT_EQUAL_VEC(expected, found, tol) {double _norm_ = std::sqrt((expected).dot(expected)); double _scale_ = _norm_ > 1.0 ? _norm_ : 1.0; if ((std::abs(((expected)[0])-((found)[0]))/_scale_ > (tol)) || (std::abs(((expected)[1])-((found)[1]))/_scale_ > (tol)) || (std::abs(((expected)[2])-((found)[2]))/_scale_ > (tol))) {std::stringstream details; details << " Expected "<<(expected)<<", found "<<(found); throwException(__FILE__, __LINE__, details.str());}};
#define ASSERT_USUALLY_TRUE(cond) {if (!(cond)) throwException(__FILE__, __LINE__, "(This test is stochastic and may occasionally fail)");};
#define ASSERT_USUALLY_EQUAL_TOL(expected, found, tol) {double _scale_ = std::abs(expected) > 1.0 ? std::abs(expected) : 1.0; if (!(std::abs((expected)-(found))/_scale_ <= (tol))) {std::stringstream details; details << "Expected "<<(expected)<<", found "<<(found)<<" (This test is stochastic and may occasionally fail)"; throwException(__FILE__, __LINE__, details.str());}};
#define ASSERT_VALID_INDEX(index, vector) {if (index < 0 || index >= (int) vector.size()) throwException(__FILE__, __LINE__, "Index out of range");};
......
openmmapi/src/MonteCarloAnisotropicBarostatImpl.cpp
View file @ 49d60231
......@@ -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) 2010 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Portions copyright (c) 2010-2013 Stanford University and the Authors. *
* Authors: Peter Eastman, Lee-Ping Wang *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
......@@ -65,7 +65,7 @@ void MonteCarloAnisotropicBarostatImpl::initialize(ContextImpl& context) {
void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context) {
if (++step < owner.getFrequency() || owner.getFrequency() == 0)
return;
if (owner.getScaleX() == 0 && owner.getScaleY() == 0 && owner.getScaleZ() == 0)
if (!owner.getScaleX() && !owner.getScaleY() && !owner.getScaleZ())
return;
step = 0;
......@@ -75,23 +75,26 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
double pressure;
// Choose which axis to modify at random.
double rnd = genrand_real2(random)*3.0;
int axis;
while (1) {
if (rnd < 1.0 && owner.getScaleX()) {
while (true) {
double rnd = genrand_real2(random)*3.0;
if (rnd < 1.0) {
if (owner.getScaleX()) {
axis = 0;
pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureX())*(AVOGADRO*1e-25);
break;
} else if (rnd < 2.0 && owner.getScaleY()) {
}
} else if (rnd < 2.0) {
if (owner.getScaleY()) {
axis = 1;
pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureY())*(AVOGADRO*1e-25);
break;
}
} else if (owner.getScaleZ()) {
axis = 2;
pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureZ())*(AVOGADRO*1e-25);
break;
}
rnd = genrand_real2(random)*3.0;
}
// Modify the periodic box size.
......@@ -101,9 +104,7 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
double volume = box[0][0]*box[1][1]*box[2][2];
double deltaVolume = volumeScale[axis]*2*(genrand_real2(random)-0.5);
double newVolume = volume+deltaVolume;
Vec3 lengthScale;
for (int i=0; i<3; i++)
lengthScale[i] = 1.0;
Vec3 lengthScale(1.0, 1.0, 1.0);
lengthScale[axis] = newVolume/volume;
kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale[0], lengthScale[1], lengthScale[2]);
context.getOwner().setPeriodicBoxVectors(box[0]*lengthScale[0], box[1]*lengthScale[1], box[2]*lengthScale[2]);
......
platforms/cuda/tests/TestCudaMonteCarloAnisotropicBarostat.cpp
View file @ 49d60231
......@@ -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-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Portions copyright (c) 2008-2013 Stanford University and the Authors. *
* Authors: Peter Eastman, Lee-Ping Wang *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
......@@ -274,42 +274,35 @@ void testRandomSeed() {
}
}
void testEinsteinCrystal() {
/*
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
/**
* 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 = 10;
const int frequency = 2;
const int equil = 10000;
const int steps = 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[] = {9.0, 10.0, 11.0};
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);
// All results.
double results[] = {0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0};
int simNum = 0;
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.
......@@ -325,6 +318,8 @@ void testEinsteinCrystal() {
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);
......@@ -332,8 +327,10 @@ void testEinsteinCrystal() {
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, frequency, (a==0||a==3), (a==1||a==3), (a==2||a==3));
system.addForce(barostat);
......@@ -352,8 +349,7 @@ void testEinsteinCrystal() {
integrator.step(frequency);
}
volume /= steps;
results[simNum] = volume;
simNum += 1;
results.push_back(volume);
}
}
for (int p = 0; p < 3; p++) {
......@@ -368,6 +364,8 @@ void testEinsteinCrystal() {
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);
......@@ -375,13 +373,15 @@ void testEinsteinCrystal() {
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.01);
LangevinIntegrator integrator(temp, 0.1, 0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
// Let it equilibrate.
......@@ -395,29 +395,22 @@ void testEinsteinCrystal() {
integrator.step(frequency);
}
volume /= steps;
results[simNum] = volume;
simNum += 1;
results.push_back(volume);
}
/*
for (int j = 0; j < 15; j++) {
printf("%.6f\n",results[j]);
}
*/
// Check to see if volumes decrease with increasing pressure.
ASSERT(results[0] > results[1]);
ASSERT(results[1] > results[2]);
ASSERT(results[3] > results[4]);
ASSERT(results[4] > results[5]);
ASSERT(results[6] > results[7]);
ASSERT(results[7] > results[8]);
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((results[0] - results[1]) > (results[3] - results[4]));
ASSERT((results[1] - results[2]) > (results[4] - results[5]));
ASSERT((results[3] - results[4]) > (results[6] - results[7]));
ASSERT((results[4] - results[5]) > (results[7] - results[8]));
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));
......
platforms/opencl/tests/TestOpenCLMonteCarloAnisotropicBarostat.cpp
View file @ 49d60231
......@@ -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-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Portions copyright (c) 2008-2013 Stanford University and the Authors. *
* Authors: Peter Eastman, Lee-Ping Wang *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
......@@ -51,7 +51,7 @@
using namespace OpenMM;
using namespace std;
static OpenCLPlatform platform;
OpenCLPlatform platform;
void testChangingBoxSize() {
System system;
......@@ -274,42 +274,35 @@ void testRandomSeed() {
}
}
void testEinsteinCrystal() {
/*
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
/**
* 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 = 10;
const int frequency = 2;
const int equil = 10000;
const int steps = 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[] = {9.0, 10.0, 11.0};
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);
// All results.
double results[] = {0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0};
int simNum = 0;
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.
......@@ -325,6 +318,8 @@ void testEinsteinCrystal() {
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);
......@@ -332,8 +327,10 @@ void testEinsteinCrystal() {
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, frequency, (a==0||a==3), (a==1||a==3), (a==2||a==3));
system.addForce(barostat);
......@@ -352,8 +349,7 @@ void testEinsteinCrystal() {
integrator.step(frequency);
}
volume /= steps;
results[simNum] = volume;
simNum += 1;
results.push_back(volume);
}
}
for (int p = 0; p < 3; p++) {
......@@ -368,6 +364,8 @@ void testEinsteinCrystal() {
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);
......@@ -375,13 +373,15 @@ void testEinsteinCrystal() {
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.01);
LangevinIntegrator integrator(temp, 0.1, 0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
// Let it equilibrate.
......@@ -395,29 +395,22 @@ void testEinsteinCrystal() {
integrator.step(frequency);
}
volume /= steps;
results[simNum] = volume;
simNum += 1;
results.push_back(volume);
}
/*
for (int j = 0; j < 15; j++) {
printf("%.6f\n",results[j]);
}
*/
// Check to see if volumes decrease with increasing pressure.
ASSERT(results[0] > results[1]);
ASSERT(results[1] > results[2]);
ASSERT(results[3] > results[4]);
ASSERT(results[4] > results[5]);
ASSERT(results[6] > results[7]);
ASSERT(results[7] > results[8]);
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((results[0] - results[1]) > (results[3] - results[4]));
ASSERT((results[1] - results[2]) > (results[4] - results[5]));
ASSERT((results[3] - results[4]) > (results[6] - results[7]));
ASSERT((results[4] - results[5]) > (results[7] - results[8]));
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));
......
platforms/reference/tests/TestReferenceMonteCarloAnisotropicBarostat.cpp
View file @ 49d60231
......@@ -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-2010 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Portions copyright (c) 2008-2013 Stanford University and the Authors. *
* Authors: Peter Eastman, Lee-Ping Wang *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
......@@ -34,6 +34,8 @@
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/MonteCarloBarostat.h"
#include "openmm/MonteCarloAnisotropicBarostat.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
......@@ -274,6 +276,151 @@ void testRandomSeed() {
}
}
/**
* 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);
ReferencePlatform platform;
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, frequency, (a==0||a==3), (a==1||a==3), (a==2||a==3));
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() {
try {
testIdealGas();
......@@ -281,6 +428,7 @@ int main() {
testIdealGasAxis(1);
testIdealGasAxis(2);
testRandomSeed();
testEinsteinCrystal();
}
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