Skip to content

GitLab

Commit 93e81862 authored by Lee-Ping's avatar Lee-Ping
Browse files

Placed MonteCarloAnisotropicBarostat into separate files

parent 3dbc4f78
Hide whitespace changes
Inline Side-by-side
Showing with 462 additions and 287 deletions
openmmapi/include/openmm/MonteCarloAnisotropicBarostat.h 0 → 100644
View file @ 93e81862
#ifndef OPENMM_MONTECARLOANISOTROPICBAROSTAT_H_
#define OPENMM_MONTECARLOANISOTROPICBAROSTAT_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 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 "Vec3.h"
#include <string>
#include "internal/windowsExport.h"
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 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 MonteCarloAnisotropicBarostat : public Force {
public:
/**
* This is the name of the parameter which stores the current pressure acting on
* the X-axis (in bar).
*/
static const std::string& PressureX() {
static const std::string key = "MonteCarloPressureX";
return key;
}
/**
* This is the name of the parameter which stores the current pressure acting on
* the Y-axis (in bar).
*/
static const std::string& PressureY() {
static const std::string key = "MonteCarloPressureY";
return key;
}
/**
* This is the name of the parameter which stores the current pressure acting on
* the Z-axis (in bar).
*/
static const std::string& PressureZ() {
static const std::string key = "MonteCarloPressureZ";
return key;
}
/**
* Create a MonteCarloAnisotropicBarostat.
*
* @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
*/
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.
*/
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.
*/
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.
*/
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.
*/
bool getScaleZ() const {
return scaleZ;
}
/**
* Get the frequency (in time steps) at which Monte Carlo pressure 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 pressure 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 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:
Vec3 defaultPressure;
double temperature;
bool scaleX, scaleY, scaleZ;
int frequency, randomNumberSeed;
double GetTemperature() const {
return temperature;
}
void SetTemperature(double temperature) {
this->temperature = temperature;
}
};
} // namespace OpenMM
#endif /*OPENMM_MONTECARLOANISOTROPICBAROSTAT_H_*/
openmmapi/include/openmm/MonteCarloBarostat.h
View file @ 93e81862
......@@ -33,7 +33,6 @@
* -------------------------------------------------------------------------- */
#include "Force.h"
#include "Vec3.h"
#include <string>
#include "internal/windowsExport.h"
......@@ -134,149 +133,6 @@ private:
}
};
/**
* 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 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 MonteCarloAnisotropicBarostat : public Force {
public:
/**
* This is the name of the parameter which stores the current pressure acting on
* the X-axis (in bar).
*/
static const std::string& PressureX() {
static const std::string key = "MonteCarloPressureX";
return key;
}
/**
* This is the name of the parameter which stores the current pressure acting on
* the Y-axis (in bar).
*/
static const std::string& PressureY() {
static const std::string key = "MonteCarloPressureY";
return key;
}
/**
* This is the name of the parameter which stores the current pressure acting on
* the Z-axis (in bar).
*/
static const std::string& PressureZ() {
static const std::string key = "MonteCarloPressureZ";
return key;
}
/**
* Create a MonteCarloAnisotropicBarostat.
*
* @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
*/
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.
*/
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.
*/
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.
*/
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.
*/
bool getScaleZ() const {
return scaleZ;
}
/**
* Get the frequency (in time steps) at which Monte Carlo pressure 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 pressure 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 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:
Vec3 defaultPressure;
double temperature;
bool scaleX, scaleY, scaleZ;
int frequency, randomNumberSeed;
double GetTemperature() const {
return temperature;
}
void SetTemperature(double temperature) {
this->temperature = temperature;
}
};
} // namespace OpenMM
#endif /*OPENMM_MONTECARLOBAROSTAT_H_*/
openmmapi/include/openmm/internal/MonteCarloAnisotropicBarostatImpl.h 0 → 100644
View file @ 93e81862
#ifndef OPENMM_MONTECARLOANISOTROPICBAROSTATIMPL_H_
#define OPENMM_MONTECARLOANISOTROPICBAROSTATIMPL_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 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/MonteCarloAnisotropicBarostat.h"
#include "openmm/Kernel.h"
#include "sfmt/SFMT.h"
#include <string>
namespace OpenMM {
/**
* This is the internal implementation of MonteCarloAnisotropicBarostat.
*/
class MonteCarloAnisotropicBarostatImpl : public ForceImpl {
public:
MonteCarloAnisotropicBarostatImpl(const MonteCarloAnisotropicBarostat& owner);
void initialize(ContextImpl& context);
const MonteCarloAnisotropicBarostat& 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 MonteCarloAnisotropicBarostat& owner;
int step, numAttempted[3], numAccepted[3];
double volumeScale[3];
OpenMM_SFMT::SFMT random;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_MONTECARLOANISOTROPICBAROSTATIMPL_H_*/
openmmapi/include/openmm/internal/MonteCarloBarostatImpl.h
View file @ 93e81862
......@@ -66,28 +66,6 @@ private:
Kernel kernel;
};
class MonteCarloAnisotropicBarostatImpl : public ForceImpl {
public:
MonteCarloAnisotropicBarostatImpl(const MonteCarloAnisotropicBarostat& owner);
void initialize(ContextImpl& context);
const MonteCarloAnisotropicBarostat& 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 MonteCarloAnisotropicBarostat& owner;
int step, numAttempted[3], numAccepted[3];
double volumeScale[3];
OpenMM_SFMT::SFMT random;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_MONTECARLOBAROSTATIMPL_H_*/
openmmapi/src/MonteCarloAnisotropicBarostat.cpp 0 → 100644
View file @ 93e81862
/* -------------------------------------------------------------------------- *
* 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 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/MonteCarloAnisotropicBarostat.h"
#include "openmm/internal/MonteCarloAnisotropicBarostatImpl.h"
#include <ctime>
using namespace OpenMM;
MonteCarloAnisotropicBarostat::MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double temperature, int frequency, bool scaleX, bool scaleY, bool scaleZ) :
defaultPressure(defaultPressure), temperature(temperature), frequency(frequency), scaleX(scaleX), scaleY(scaleY), scaleZ(scaleZ) {
setRandomNumberSeed((int) time(NULL));
}
ForceImpl* MonteCarloAnisotropicBarostat::createImpl() const {
return new MonteCarloAnisotropicBarostatImpl(*this);
}
openmmapi/src/MonteCarloAnisotropicBarostatImpl.cpp 0 → 100644
View file @ 93e81862
/* -------------------------------------------------------------------------- *
* 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 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/MonteCarloAnisotropicBarostatImpl.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/Context.h"
#include "openmm/kernels.h"
#include <cmath>
#include <vector>
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))
MonteCarloAnisotropicBarostatImpl::MonteCarloAnisotropicBarostatImpl(const MonteCarloAnisotropicBarostat& owner) : owner(owner), step(0) {
}
void MonteCarloAnisotropicBarostatImpl::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 MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context) {
if (++step < owner.getFrequency() || owner.getFrequency() == 0)
return;
if (owner.getScaleX() == 0 && owner.getScaleY() == 0 && owner.getScaleZ() == 0)
return;
step = 0;
// Compute the current potential energy.
double initialEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
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()) {
axis = 0;
pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureX())*(AVOGADRO*1e-25);
break;
} else if (rnd < 2.0 && 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;
}
}
// 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;
for (int i=0; i<3; i++)
lengthScale[i] = 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]);
// 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 - 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> MonteCarloAnisotropicBarostatImpl::getDefaultParameters() {
std::map<std::string, double> parameters;
parameters[MonteCarloAnisotropicBarostat::PressureX()] = getOwner().getDefaultPressure()[0];
parameters[MonteCarloAnisotropicBarostat::PressureY()] = getOwner().getDefaultPressure()[1];
parameters[MonteCarloAnisotropicBarostat::PressureZ()] = getOwner().getDefaultPressure()[2];
return parameters;
}
std::vector<std::string> MonteCarloAnisotropicBarostatImpl::getKernelNames() {
std::vector<std::string> names;
names.push_back(ApplyMonteCarloBarostatKernel::Name());
return names;
}
openmmapi/src/MonteCarloBarostat.cpp
View file @ 93e81862
......@@ -43,12 +43,3 @@ MonteCarloBarostat::MonteCarloBarostat(double defaultPressure, double temperatur
ForceImpl* MonteCarloBarostat::createImpl() const {
return new MonteCarloBarostatImpl(*this);
}
MonteCarloAnisotropicBarostat::MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double temperature, int frequency, bool scaleX, bool scaleY, bool scaleZ) :
defaultPressure(defaultPressure), temperature(temperature), frequency(frequency), scaleX(scaleX), scaleY(scaleY), scaleZ(scaleZ) {
setRandomNumberSeed((int) time(NULL));
}
ForceImpl* MonteCarloAnisotropicBarostat::createImpl() const {
return new MonteCarloAnisotropicBarostatImpl(*this);
}
openmmapi/src/MonteCarloBarostatImpl.cpp
View file @ 93e81862
......@@ -64,22 +64,22 @@ void MonteCarloBarostatImpl::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();
// 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*2*(genrand_real2(random)-0.5);
double newVolume = volume+deltaVolume;
double lengthScale = std::pow(newVolume/volume, 1.0/3.0);
kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale, lengthScale, lengthScale);
kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale);
context.getOwner().setPeriodicBoxVectors(box[0]*lengthScale, box[1]*lengthScale, box[2]*lengthScale);
// Compute the energy of the modified system.
double finalEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
......@@ -88,7 +88,7 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context) {
double w = finalEnergy-initialEnergy + pressure*deltaVolume - 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;
......@@ -122,108 +122,3 @@ std::vector<std::string> MonteCarloBarostatImpl::getKernelNames() {
return names;
}
MonteCarloAnisotropicBarostatImpl::MonteCarloAnisotropicBarostatImpl(const MonteCarloAnisotropicBarostat& owner) : owner(owner), step(0) {
}
void MonteCarloAnisotropicBarostatImpl::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 MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context) {
if (++step < owner.getFrequency() || owner.getFrequency() == 0)
return;
if (owner.getScaleX() == 0 && owner.getScaleY() == 0 && owner.getScaleZ() == 0)
return;
step = 0;
// Compute the current potential energy.
double initialEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
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()) {
axis = 0;
pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureX())*(AVOGADRO*1e-25);
break;
} else if (rnd < 2.0 && 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;
}
}
// 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;
for (int i=0; i<3; i++)
lengthScale[i] = 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]);
// 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 - 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> MonteCarloAnisotropicBarostatImpl::getDefaultParameters() {
std::map<std::string, double> parameters;
parameters[MonteCarloAnisotropicBarostat::PressureX()] = getOwner().getDefaultPressure()[0];
parameters[MonteCarloAnisotropicBarostat::PressureY()] = getOwner().getDefaultPressure()[1];
parameters[MonteCarloAnisotropicBarostat::PressureZ()] = getOwner().getDefaultPressure()[2];
return parameters;
}
std::vector<std::string> MonteCarloAnisotropicBarostatImpl::getKernelNames() {
std::vector<std::string> names;
names.push_back(ApplyMonteCarloBarostatKernel::Name());
return names;
}
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