Commit 5c7321ca authored by peastman's avatar peastman
Browse files

Created RPMDMonteCarloBarostat

parent 1c271535
#ifndef OPENMM_RPMDMONTECARLOBAROSTAT_H_
#define OPENMM_RPMDMONTECARLOBAROSTAT_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-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/Force.h"
#include <string>
#include "internal/windowsExportRpmd.h"
namespace OpenMM {
/**
* This class is very similar to MonteCarloBarostat, but it is specifically designed for use
* with RPMDIntegrator. For each trial move, it scales all copies of the system by the same
* amount, then accepts or rejects the move based on the change to the total energy of the
* ring polymer (as returned by the integrator's getTotalEnergy() method).
*/
class OPENMM_EXPORT_RPMD RPMDMonteCarloBarostat : public Force {
public:
/**
* 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 = "RPMDMonteCarloPressure";
return key;
}
/**
* Create a MonteCarloBarostat.
*
* @param defaultPressure the default pressure acting on the system (in bar)
* @param frequency the frequency at which Monte Carlo pressure changes should be attempted (in time steps)
*/
RPMDMonteCarloBarostat(double defaultPressure, 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 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 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.
*
* If seed is set to 0 (which is the default value assigned), a unique seed is chosen when a Context
* is created from this Force. This is done to ensure that each Context receives unique random seeds
* without you needing to set them explicitly.
*/
void setRandomNumberSeed(int seed) {
randomNumberSeed = seed;
}
/**
* Returns whether or not this force makes use of periodic boundary
* conditions.
*
* @returns true if force uses PBC and false otherwise
*/
bool usesPeriodicBoundaryConditions() const {
return true;
}
protected:
ForceImpl* createImpl() const;
private:
double defaultPressure;
int frequency, randomNumberSeed;
};
} // namespace OpenMM
#endif /*OPENMM_RPMDMONTECARLOBAROSTAT_H_*/
#ifndef OPENMM_RPMDUPDATER_H_
#define OPENMM_RPMDUPDATER_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) 2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/ForceImpl.h"
#include "internal/windowsExportRpmd.h"
namespace OpenMM {
/**
* This ForceImpl subclass is for forces specifically designed to work with RPMDIntegrator.
* It adds an updateRPMDState() method that is invoked at the start of each time step
* (in contrast to updateContextState(), which gets invoked many times per time step, once
* for each copy).
*/
class OPENMM_EXPORT_RPMD RPMDUpdater : public ForceImpl {
public:
virtual void updateRPMDState(ContextImpl& context) = 0;
};
} // namespace OpenMM
#endif /*OPENMM_RPMDUPDATER_H_*/
#ifndef OPENMM_RPMDMONTECARLOBAROSTATIMPL_H_
#define OPENMM_RPMDMONTECARLOBAROSTATIMPL_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-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/RPMDMonteCarloBarostat.h"
#include "openmm/RPMDUpdater.h"
#include "openmm/Kernel.h"
#include "openmm/Vec3.h"
#include "sfmt/SFMT.h"
#include <string>
#include <vector>
namespace OpenMM {
/**
* This is the internal implementation of RPMDMonteCarloBarostat.
*/
class RPMDMonteCarloBarostatImpl : public RPMDUpdater {
public:
RPMDMonteCarloBarostatImpl(const RPMDMonteCarloBarostat& owner);
void initialize(ContextImpl& context);
const RPMDMonteCarloBarostat& getOwner() const {
return owner;
}
void updateRPMDState(ContextImpl& context);
void updateContextState(ContextImpl& context) {
// This is unused, since the updating is done in updateRPMDState().
}
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 RPMDMonteCarloBarostat& owner;
int step, numAttempted, numAccepted;
double volumeScale;
OpenMM_SFMT::SFMT random;
std::vector<std::vector<Vec3> > savedPositions;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_RPMDMONTECARLOBAROSTATIMPL_H_*/
......@@ -34,6 +34,7 @@
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/RpmdKernels.h"
#include "openmm/RPMDUpdater.h"
#include "SimTKOpenMMRealType.h"
#include <cmath>
#include <string>
......@@ -190,6 +191,12 @@ void RPMDIntegrator::step(int steps) {
context->getOwner().setPositions(p);
isFirstStep = false;
}
vector<ForceImpl*>& forceImpls = context->getForceImpls();
for (int i = 0; i < (int) forceImpls.size(); i++) {
RPMDUpdater* updater = dynamic_cast<RPMDUpdater*>(forceImpls[i]);
if (updater != NULL)
updater->updateRPMDState(*context);
}
for (int i = 0; i < steps; ++i) {
kernel.getAs<IntegrateRPMDStepKernel>().execute(*context, *this, forcesAreValid);
forcesAreValid = true;
......
/* -------------------------------------------------------------------------- *
* 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-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/RPMDMonteCarloBarostat.h"
#include "openmm/internal/RPMDMonteCarloBarostatImpl.h"
using namespace OpenMM;
RPMDMonteCarloBarostat::RPMDMonteCarloBarostat(double defaultPressure, int frequency) :
defaultPressure(defaultPressure), frequency(frequency) {
setRandomNumberSeed(0);
}
ForceImpl* RPMDMonteCarloBarostat::createImpl() const {
return new RPMDMonteCarloBarostatImpl(*this);
}
/* -------------------------------------------------------------------------- *
* 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-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/RPMDMonteCarloBarostatImpl.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/OSRngSeed.h"
#include "openmm/Context.h"
#include "openmm/kernels.h"
#include "openmm/OpenMMException.h"
#include "openmm/RPMDIntegrator.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))
RPMDMonteCarloBarostatImpl::RPMDMonteCarloBarostatImpl(const RPMDMonteCarloBarostat& owner) : owner(owner), step(0) {
}
void RPMDMonteCarloBarostatImpl::initialize(ContextImpl& context) {
RPMDIntegrator* integrator = dynamic_cast<RPMDIntegrator*>(&context.getIntegrator());
if (integrator == NULL)
throw OpenMMException("RPMDMonteCarloBarostat must be used with an RPMDIntegrator");;
if (!integrator->getApplyThermostat())
throw OpenMMException("RPMDMonteCarloBarostat requires the integrator's thermostat to be enabled");;
kernel = context.getPlatform().createKernel(ApplyMonteCarloBarostatKernel::Name(), context);
kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner);
savedPositions.resize(integrator->getNumCopies());
Vec3 box[3];
context.getPeriodicBoxVectors(box[0], box[1], box[2]);
double volume = box[0][0]*box[1][1]*box[2][2];
volumeScale = 0.01*volume;
numAttempted = 0;
numAccepted = 0;
int randSeed = owner.getRandomNumberSeed();
// A random seed of 0 means use a unique one
if (randSeed == 0) randSeed = osrngseed();
init_gen_rand(randSeed, random);
}
void RPMDMonteCarloBarostatImpl::updateRPMDState(ContextImpl& context) {
if (++step < owner.getFrequency() || owner.getFrequency() == 0)
return;
step = 0;
// Compute the current potential energy.
RPMDIntegrator& integrator = dynamic_cast<RPMDIntegrator&>(context.getIntegrator());
// Record the initial positions and energy
double initialEnergy = 0;
int numCopies = integrator.getNumCopies();
for (int i = 0; i < numCopies; i++) {
State state = integrator.getState(i, State::Positions | State::Energy);
savedPositions[i] = state.getPositions();
initialEnergy += state.getPotentialEnergy();
}
// Compute the centroid.
int numParticles = context.getSystem().getNumParticles();
vector<Vec3> centroid(numParticles, Vec3());
for (int i = 0; i < numParticles; i++) {
for (int j = 0; j < numCopies; j++)
centroid[i] += savedPositions[j][i];
centroid[i] *= 1.0/numCopies;
}
// Modify the periodic box size and scale the coordinates of the centroid.
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);
context.setPositions(centroid);
kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale, lengthScale, lengthScale);
context.getOwner().setPeriodicBoxVectors(box[0]*lengthScale, box[1]*lengthScale, box[2]*lengthScale);
State scaledState = context.getOwner().getState(State::Positions);
// Now apply the same offset to all the copies.
vector<Vec3> delta(numParticles);
for (int i = 0; i < numParticles; i++)
delta[i] = scaledState.getPositions()[i]-centroid[i];
double finalEnergy = 0;
vector<Vec3> positions(numParticles);
for (int copy = 0; copy < numCopies; copy++) {
for (int i = 0; i < numParticles; i++)
positions[i] = savedPositions[copy][i]+delta[i];
integrator.setPositions(copy, positions);
finalEnergy += integrator.getState(copy, State::Energy).getPotentialEnergy();
}
// Compute the energy of the modified system.
double pressure = context.getParameter(RPMDMonteCarloBarostat::Pressure())*(AVOGADRO*1e-25);
double kT = BOLTZ*integrator.getTemperature();
double w = (finalEnergy-initialEnergy)/numCopies + pressure*deltaVolume - context.getMolecules().size()*kT*std::log(newVolume/volume);
if (w > 0 && genrand_real2(random) > std::exp(-w/kT)) {
// Reject the step.
for (int copy = 0; copy < numCopies; copy++)
integrator.setPositions(copy, savedPositions[copy]);
context.getOwner().setPeriodicBoxVectors(box[0], box[1], box[2]);
volume = newVolume;
}
else
numAccepted++;
numAttempted++;
if (numAttempted >= 10) {
if (numAccepted < 0.25*numAttempted) {
volumeScale /= 1.1;
numAttempted = 0;
numAccepted = 0;
}
else if (numAccepted > 0.75*numAttempted) {
volumeScale = std::min(volumeScale*1.1, volume*0.3);
numAttempted = 0;
numAccepted = 0;
}
}
}
std::map<std::string, double> RPMDMonteCarloBarostatImpl::getDefaultParameters() {
std::map<std::string, double> parameters;
parameters[RPMDMonteCarloBarostat::Pressure()] = getOwner().getDefaultPressure();
return parameters;
}
std::vector<std::string> RPMDMonteCarloBarostatImpl::getKernelNames() {
std::vector<std::string> names;
names.push_back(ApplyMonteCarloBarostatKernel::Name());
return names;
}
......@@ -48,6 +48,7 @@ namespace std {
#include "OpenMM.h"
#include "OpenMMAmoeba.h"
#include "openmm/RPMDIntegrator.h"
#include "openmm/RPMDMonteCarloBarostat.h"
#include "OpenMMDrude.h"
#include "openmm/serialization/SerializationNode.h"
#include "openmm/serialization/SerializationProxy.h"
......
......@@ -621,6 +621,7 @@ EXCLUDE_PATTERNS = */tests/* \
*amoebaKernels.h \
*DrudeKernels.h \
*RpmdKernels.h \
*RPMDUpdater.h \
*OpenMMFortranModule.f90 \
*OpenMMCWrapper.h
......
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