Skip to content

GitLab

Commit 7d1e3fc4 authored by peastman's avatar peastman
Browse files

Created MonteCarloMembraneBarostat

parent 4c5f80ce
Show whitespace changes
Inline Side-by-side
Showing with 706 additions and 0 deletions
openmmapi/include/openmm/MonteCarloMembraneBarostat.h 0 → 100644
View file @ 7d1e3fc4
#ifndef OPENMM_MONTECARLOMEMBRANEBAROSTAT_H_
#define OPENMM_MONTECARLOMEMBRANEBAROSTAT_H_
/* -------------------------------------------------------------------------- *
* 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) 2010-2014 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 "Force.h"
#include <string>
#include "internal/windowsExport.h"
namespace OpenMM {
/**
* This is a Monte Carlo barostat designed specifically for membrane simulations. It assumes the
* membrane lies in the XY plane. The Monte Carlo acceptance criterion includes a term to model
* isotropic pressure, which depends on the volume of the periodic box, and a second term to model
* surface tension, which depends on the cross sectional area of the box in the XY plane. Note
* that pressure and surface tension are defined with opposite senses: a larger pressure tends to
* make the box smaller, but a larger surface tension tends to make the box larger.
*
* There are options for configuring exactly how the various box dimensions are allowed to change:
*
* <ul>
* <li>The X and Y axes may be treated isotropically, in which case they always scale by the same
* amount and remain in proportion to each other; or they may be treated anisotropically, in which
* case they can vary independently of each other.</li>
* <li>The Z axis can be allowed to vary independently of the other axes; or held fixed; or constrained
* to vary in inverse proportion to the other two axes, so that the total box volume remains fixed.</li>
* </ul>
*
* 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
* moves). It does not actually perform temperature regulation, however. You must use another
* mechanism along with it to maintain the temperature, such as LangevinIntegrator or AndersenThermostat.
*/
class OPENMM_EXPORT MonteCarloMembraneBarostat : public Force {
public:
/**
* This is an enumeration of the different behaviors for the X and Y axes.
*/
enum XYMode {
/**
* The X and Y axes are always scaled by the same amount, so the ratio of their lengths remains constant.
*/
XYIsotropic = 0,
/**
* The X and Y axes are allowed to vary independently of each other.
*/
XYAnisotropic = 1
};
/**
* This is an enumeration of the different behaviors for Z axis.
*/
enum ZMode {
/**
* The Z axis is allowed to vary freely, independent of the other two axes.
*/
ZFree = 0,
/**
* The Z axis is held fixed and does not change.
*/
ZFixed = 1,
/**
* The Z axis is always scaled in inverse proportion to the other two axes so the box volume remains
* fixed. Note that in this mode pressure has no effect on the system, only surface tension.
*/
ConstantVolume = 2
};
/**
* This is the name of the parameter which stores the current pressure acting on
* the system (in bar).
*/
static const std::string& Pressure() {
static const std::string key = "MembraneMonteCarloPressure";
return key;
}
/**
* This is the name of the parameter which stores the current surface tension acting on
* the system (in bar*nm).
*/
static const std::string& SurfaceTension() {
static const std::string key = "MembraneMonteCarloSurfaceTension";
return key;
}
/**
* Create a MonteCarloMembraneBarostat.
*
* @param defaultPressure the default pressure acting on the system (in bar)
* @param defaultSurfaceTension the default surface tension acting on the system (in bar*nm)
* @param temperature the temperature at which the system is being maintained (in Kelvin)
* @param xymode the mode specifying the behavior of the X and Y axes
* @param zmode the mode specifying the behavior of the Z axis
* @param frequency the frequency at which Monte Carlo volume changes should be attempted (in time steps)
*/
MonteCarloMembraneBarostat(double defaultPressure, double defaultSurfaceTension, double temperature, XYMode xymode, ZMode zmode, int frequency = 25);
/**
* Get the default pressure acting on the system (in bar).
*
* @return the default pressure acting on the system, measured in bar.
*/
double getDefaultPressure() const {
return defaultPressure;
}
/**
* Set the default pressure acting on the system. This will affect any new Contexts you create,
* but not ones that already exist.
*
* @param pressure the default pressure acting on the system, measured in bar.
*/
void setDefaultPressure(double pressure) {
defaultPressure = pressure;
}
/**
* Get the default surface tension acting on the system (in bar*nm).
*
* @return the default surface tension acting on the system, measured in bar*nm.
*/
double getDefaultSurfaceTension() const {
return defaultSurfaceTension;
}
/**
* Set the default surface tension acting on the system. This will affect any new Contexts you create,
* but not ones that already exist.
*
* @param surfaceTension the default surface tension acting on the system, measured in bar.
*/
void setDefaultSurfaceTension(double surfaceTension) {
defaultSurfaceTension = surfaceTension;
}
/**
* Get the frequency (in time steps) at which Monte Carlo volume changes should be attempted. If this is set to
* 0, the barostat is disabled.
*/
int getFrequency() const {
return frequency;
}
/**
* Set the frequency (in time steps) at which Monte Carlo volume changes should be attempted. If this is set to
* 0, the barostat is disabled.
*/
void setFrequency(int freq) {
frequency = freq;
}
/**
* Get the temperature at which the system is being maintained, measured in Kelvin.
*/
double getTemperature() const {
return temperature;
}
/**
* Set the temperature at which the system is being maintained.
*
* @param temp the system temperature, measured in Kelvin.
*/
void setTemperature(double temp) {
temperature = temp;
}
/**
* Get the mode specifying the behavior of the X and Y axes.
*/
XYMode getXYMode() const {
return xymode;
}
/**
* Set the mode specifying the behavior of the X and Y axes.
*/
void setXYMode(XYMode mode) {
xymode = mode;
}
/**
* Get the mode specifying the behavior of the Z axis.
*/
ZMode getZMode() const {
return zmode;
}
/**
* Set the mode specifying the behavior of the Z axis.
*/
void setZMode(ZMode mode) {
zmode = mode;
}
/**
* Get the random number seed. See setRandomNumberSeed() for details.
*/
int getRandomNumberSeed() const {
return randomNumberSeed;
}
/**
* Set the random number seed. It is guaranteed that if two simulations are run
* with different random number seeds, the sequence of Monte Carlo steps will be different. On
* the other hand, no guarantees are made about the behavior of simulations that use the same seed.
* In particular, Platforms are permitted to use non-deterministic algorithms which produce different
* results on successive runs, even if those runs were initialized identically.
*/
void setRandomNumberSeed(int seed) {
randomNumberSeed = seed;
}
protected:
ForceImpl* createImpl() const;
private:
double defaultPressure, defaultSurfaceTension, temperature;
XYMode xymode;
ZMode zmode;
int frequency, randomNumberSeed;
};
} // namespace OpenMM
#endif /*OPENMM_MONTECARLOMEMBRANEBAROSTAT_H_*/
openmmapi/include/openmm/internal/MonteCarloMembraneBarostatImpl.h 0 → 100644
View file @ 7d1e3fc4
#ifndef OPENMM_MONTECARLOMEMBRANEBAROSTATIMPL_H_
#define OPENMM_MONTECARLOMEMBRANEBAROSTATIMPL_H_
/* -------------------------------------------------------------------------- *
* 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) 2010-2014 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 "ForceImpl.h"
#include "openmm/MonteCarloMembraneBarostat.h"
#include "openmm/Kernel.h"
#include "sfmt/SFMT.h"
#include <string>
namespace OpenMM {
/**
* This is the internal implementation of MonteCarloMembraneBarostat.
*/
class MonteCarloMembraneBarostatImpl : public ForceImpl {
public:
MonteCarloMembraneBarostatImpl(const MonteCarloMembraneBarostat& owner);
void initialize(ContextImpl& context);
const MonteCarloMembraneBarostat& getOwner() const {
return owner;
}
void updateContextState(ContextImpl& context);
double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
// This force doesn't apply forces to particles.
return 0.0;
}
std::map<std::string, double> getDefaultParameters();
std::vector<std::string> getKernelNames();
private:
const MonteCarloMembraneBarostat& owner;
int step, numAttempted[3], numAccepted[3];
double volumeScale[3];
OpenMM_SFMT::SFMT random;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_MONTECARLOMEMBRANEBAROSTATIMPL_H_*/
openmmapi/src/MonteCarloMembraneBarostat.cpp 0 → 100644
View file @ 7d1e3fc4
/* -------------------------------------------------------------------------- *
* 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) 2010-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/MonteCarloMembraneBarostat.h"
#include "openmm/internal/MonteCarloMembraneBarostatImpl.h"
#include "openmm/internal/OSRngSeed.h"
using namespace OpenMM;
MonteCarloMembraneBarostat::MonteCarloMembraneBarostat(double defaultPressure, double defaultSurfaceTension, double temperature, XYMode xymode, ZMode zmode, int frequency) :
defaultPressure(defaultPressure), defaultSurfaceTension(defaultSurfaceTension), temperature(temperature),
xymode(xymode), zmode(zmode), frequency(frequency) {
setRandomNumberSeed(osrngseed());
}
ForceImpl* MonteCarloMembraneBarostat::createImpl() const {
return new MonteCarloMembraneBarostatImpl(*this);
}
openmmapi/src/MonteCarloMembraneBarostatImpl.cpp 0 → 100644
View file @ 7d1e3fc4
/* -------------------------------------------------------------------------- *
* 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) 2010-2014 Stanford University and the Authors. *
* Authors: Peter Eastman, Lee-Ping Wang *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/MonteCarloMembraneBarostatImpl.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/Context.h"
#include "openmm/kernels.h"
#include <cmath>
#include <vector>
#include <algorithm>
using namespace OpenMM;
using namespace OpenMM_SFMT;
using std::vector;
const float BOLTZMANN = 1.380658e-23f; // (J/K)
const float AVOGADRO = 6.0221367e23f;
const float RGAS = BOLTZMANN*AVOGADRO; // (J/(mol K))
const float BOLTZ = RGAS/1000; // (kJ/(mol K))
MonteCarloMembraneBarostatImpl::MonteCarloMembraneBarostatImpl(const MonteCarloMembraneBarostat& owner) : owner(owner), step(0) {
}
void MonteCarloMembraneBarostatImpl::initialize(ContextImpl& context) {
kernel = context.getPlatform().createKernel(ApplyMonteCarloBarostatKernel::Name(), context);
kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner);
Vec3 box[3];
context.getPeriodicBoxVectors(box[0], box[1], box[2]);
double volume = box[0][0]*box[1][1]*box[2][2];
for (int i=0; i<3; i++) {
volumeScale[i] = 0.01*volume;
numAttempted[i] = 0;
numAccepted[i] = 0;
}
init_gen_rand(owner.getRandomNumberSeed(), random);
}
void MonteCarloMembraneBarostatImpl::updateContextState(ContextImpl& context) {
if (++step < owner.getFrequency() || owner.getFrequency() == 0)
return;
step = 0;
// Compute the current potential energy.
double initialEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
double pressure = context.getParameter(MonteCarloMembraneBarostat::Pressure())*(AVOGADRO*1e-25);
double tension = context.getParameter(MonteCarloMembraneBarostat::SurfaceTension())*(AVOGADRO*1e-25);
// Choose which axis to modify at random.
int axis;
while (true) {
double rnd = genrand_real2(random)*3.0;
if (rnd < 1.0) {
axis = 0;
break;
} else if (rnd < 2.0) {
axis = (owner.getXYMode() == MonteCarloMembraneBarostat::XYIsotropic ? 0 : 1);
break;
} else if (owner.getZMode() == MonteCarloMembraneBarostat::ZFree) {
axis = 2;
break;
}
}
// Modify the periodic box size.
Vec3 box[3];
context.getPeriodicBoxVectors(box[0], box[1], box[2]);
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(1.0, 1.0, 1.0);
if ((axis == 0 || axis == 1) && owner.getXYMode() == MonteCarloMembraneBarostat::XYIsotropic)
lengthScale[0] = lengthScale[1] = sqrt(newVolume/volume);
else
lengthScale[axis] = newVolume/volume;
if (owner.getZMode() == MonteCarloMembraneBarostat::ConstantVolume) {
lengthScale[2] = 1.0/(lengthScale[0]*lengthScale[1]);
newVolume = volume;
deltaVolume = 0;
}
double deltaArea = box[0][0]*lengthScale[0]*box[1][1]*lengthScale[1] - box[0][0]*box[1][1];
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]);
// Compute the energy of the modified system.
double finalEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
double kT = BOLTZ*owner.getTemperature();
double w = finalEnergy-initialEnergy + pressure*deltaVolume - tension*deltaArea - context.getMolecules().size()*kT*std::log(newVolume/volume);
if (w > 0 && genrand_real2(random) > std::exp(-w/kT)) {
// Reject the step.
kernel.getAs<ApplyMonteCarloBarostatKernel>().restoreCoordinates(context);
context.getOwner().setPeriodicBoxVectors(box[0], box[1], box[2]);
volume = newVolume;
}
else
numAccepted[axis]++;
numAttempted[axis]++;
if (numAttempted[axis] >= 10) {
if (numAccepted[axis] < 0.25*numAttempted[axis]) {
volumeScale[axis] /= 1.1;
numAttempted[axis] = 0;
numAccepted[axis] = 0;
}
else if (numAccepted[axis] > 0.75*numAttempted[axis]) {
volumeScale[axis] = std::min(volumeScale[axis]*1.1, volume*0.3);
numAttempted[axis] = 0;
numAccepted[axis] = 0;
}
}
}
std::map<std::string, double> MonteCarloMembraneBarostatImpl::getDefaultParameters() {
std::map<std::string, double> parameters;
parameters[MonteCarloMembraneBarostat::Pressure()] = getOwner().getDefaultPressure();
parameters[MonteCarloMembraneBarostat::SurfaceTension()] = getOwner().getDefaultSurfaceTension();
return parameters;
}
std::vector<std::string> MonteCarloMembraneBarostatImpl::getKernelNames() {
std::vector<std::string> names;
names.push_back(ApplyMonteCarloBarostatKernel::Name());
return names;
}
platforms/reference/tests/TestReferenceMonteCarloMembraneBarostat.cpp 0 → 100644
View file @ 7d1e3fc4
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2014 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. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of MonteCarloMembraneBarostat.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/MonteCarloMembraneBarostat.h"
#include "openmm/Context.h"
#include "ReferencePlatform.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;
void testIdealGas(MonteCarloMembraneBarostat::XYMode xymode, MonteCarloMembraneBarostat::ZMode zmode) {
const int numParticles = 64;
const int frequency = 1;
const int steps = 5000;
const double pressure = 1.5;
const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
const double tension = (zmode == MonteCarloMembraneBarostat::ZFixed ? 0.2 : 0.0);
const double tensionInMD = tension*(AVOGADRO*1e-25); // surface tension in kJ/mol/nm^2
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.
ReferencePlatform platform;
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));
}
MonteCarloMembraneBarostat* barostat = new MonteCarloMembraneBarostat(pressure, tension, temp[0], xymode, zmode, 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(1000);
// Now run it for a while and see if the volume is correct.
double volume = 0.0, zsize = 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];
zsize += box[2][2];
if (xymode == MonteCarloMembraneBarostat::XYIsotropic)
ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
if (zmode == MonteCarloMembraneBarostat::ZFixed)
ASSERT_EQUAL_TOL(2*initialLength, box[2][2], 1e-5);
if (zmode == MonteCarloMembraneBarostat::ConstantVolume)
ASSERT_EQUAL_TOL(initialVolume, box[0][0]*box[1][1]*box[2][2], 1e-5);
integrator.step(frequency);
}
volume /= steps;
zsize /= steps;
if (zmode != MonteCarloMembraneBarostat::ConstantVolume) {
double effectivePressure = pressureInMD-tensionInMD/zsize;
double expected = (numParticles+1)*BOLTZ*temp[i]/effectivePressure;
ASSERT_USUALLY_EQUAL_TOL(expected, volume, 0.05);
}
}
}
void testRandomSeed() {
const int numParticles = 8;
const double temp = 100.0;
const double pressure = 1.5;
const double tension = 0.3;
ReferencePlatform platform;
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);
MonteCarloMembraneBarostat* barostat = new MonteCarloMembraneBarostat(pressure, tension, temp, MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ZFree, 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]);
}
}
}
int main() {
try {
testIdealGas(MonteCarloMembraneBarostat::XYIsotropic, MonteCarloMembraneBarostat::ZFree);
testIdealGas(MonteCarloMembraneBarostat::XYIsotropic, MonteCarloMembraneBarostat::ZFixed);
testIdealGas(MonteCarloMembraneBarostat::XYIsotropic, MonteCarloMembraneBarostat::ConstantVolume);
testIdealGas(MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ZFree);
testIdealGas(MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ZFixed);
testIdealGas(MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ConstantVolume);
testRandomSeed();
}
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