Skip to content

GitLab

Commit dfba24ea authored by Peter Eastman's avatar Peter Eastman
Browse files

Created CustomIntegrator, including reference implementation

parent e9c5ab17
Show whitespace changes
Inline Side-by-side
Showing with 1461 additions and 2 deletions
olla/include/openmm/kernels.h
View file @ dfba24ea
......@@ -41,6 +41,7 @@
#include "openmm/CustomExternalForce.h"
#include "openmm/CustomGBForce.h"
#include "openmm/CustomHbondForce.h"
#include "openmm/CustomIntegrator.h"
#include "openmm/CustomNonbondedForce.h"
#include "openmm/CustomTorsionForce.h"
#include "openmm/GBSAOBCForce.h"
......@@ -143,7 +144,7 @@ public:
/**
* Set the positions of all particles.
*
* @param positions a vector containg the particle positions
* @param positions a vector containing the particle positions
*/
virtual void setPositions(ContextImpl& context, const std::vector<Vec3>& positions) = 0;
/**
......@@ -155,7 +156,7 @@ public:
/**
* Set the velocities of all particles.
*
* @param velocities a vector containg the particle velocities
* @param velocities a vector containing the particle velocities
*/
virtual void setVelocities(ContextImpl& context, const std::vector<Vec3>& velocities) = 0;
/**
......@@ -782,6 +783,66 @@ public:
virtual void execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) = 0;
};
/**
* This kernel is invoked by CustomIntegrator to take one time step.
*/
class IntegrateCustomStepKernel : public KernelImpl {
public:
static std::string Name() {
return "IntegrateCustomStep";
}
IntegrateCustomStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param integrator the CustomIntegrator this kernel will be used for
*/
virtual void initialize(const System& system, const CustomIntegrator& integrator) = 0;
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the CustomIntegrator this kernel is being used for
* @param forcesAreValid if the context has been modified since the last time step, this will be
* false to show that cached forces are invalid and must be recalculated.
* On exit, this should specify whether the cached forces are valid at the
* end of the step.
*/
virtual void execute(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid) = 0;
/**
* Get the values of all global variables.
*
* @param context the context in which to execute this kernel
* @param values on exit, this contains the values
*/
virtual void getGlobalVariables(ContextImpl& context, std::vector<double>& values) const = 0;
/**
* Set the values of all global variables.
*
* @param context the context in which to execute this kernel
* @param values a vector containing the values
*/
virtual void setGlobalVariables(ContextImpl& context, const std::vector<double>& values) = 0;
/**
* Get the values of a per-DOF variable.
*
* @param context the context in which to execute this kernel
* @param variable the index of the variable to get
* @param values on exit, this contains the values
*/
virtual void getPerDofVariable(ContextImpl& context, int variable, std::vector<Vec3>& values) const = 0;
/**
* Set the values of a per-DOF variable.
*
* @param context the context in which to execute this kernel
* @param variable the index of the variable to get
* @param values a vector containing the values
*/
virtual void setPerDofVariable(ContextImpl& context, int variable, const std::vector<Vec3>& values) = 0;
};
/**
* This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
*/
......
openmmapi/include/openmm/CustomIntegrator.h 0 → 100644
View file @ dfba24ea
#ifndef OPENMM_CUSTOMINTEGRATOR_H_
#define OPENMM_CUSTOMINTEGRATOR_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) 2011 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 "Integrator.h"
#include "Vec3.h"
#include "openmm/Kernel.h"
#include "internal/windowsExport.h"
#include <string>
#include <vector>
namespace OpenMM {
/**
* This is an Integrator that can be used to implemented arbitrary, user defined
* integration algorithms. It is flexible enough to support a wide range of
* methods including both deterministic and stochastic integrators, Metropolized
* integrators, and integrators that must integrate additional quantities along
* with the particle positions and momenta.
*
* To create an integration algorithm, you first define a set of variables the
* integrator will compute. Variables come in two types: <i>global</i> variables
* have a single value, while <i>per-DOF</i> variables have a value for every
* degree of freedom (x, y, or z coordinate of a particle). You can define as
* many variables as you want of each type. The value of any variable can be
* computed by the integration algorithm, or set directly by calling a method on
* the CustomIntegrator. All variables are persistent between integration
* steps; once a value is set, it keeps that value until it is changed by the
* user or recomputed in a later integration step.
*
* Next, you define the algorithm as a series of computations. To execute a
* time step, the integrator performs the list of computations in order. Each
* computation updates the value of one global or per-DOF value. There are
* several types of computations that can be done:
*
* <ul>
* <li>Global: You provide a mathematical expression involving only global
* variables. It is evaluated and stored into a global variable.</li>
* <li>Per-DOF: You provide a mathematical expression involving both global and
* per-DOF variables. It is evaluated once for every degree of freedom, and
* the values are stored into a per-DOF variable.</li>
* <li>Sum: You provide a mathematical expression involving both global and
* per-DOF variables. It is evaluated once for every degree of freedom. All
* of those values are then added together, and the sum is stored into a global
* variable.</li>
* <li>Constrain Positions: The particle positions are updated so that all
* distance constraints are satisfied.</li>
* <li>Constrain Velocities: The particle velocities are updated so the net
* velocity along any constrained distance is 0.</li>
* </ul>
*
* In addition to the variables you define by calling addGlobalVariable() and
* addPerDofVariable(), the integrator provides the following pre-defined
* variables:
*
* <ul>
* <li>dt: (global) This is the step size being used by the integrator.</li>
* <li>energy: (global, read-only) This is the current potential energy of the
* system.</li>
* <li>x: (per-DOF) This is the current value of the degree of freedom (the x,
* y, or z coordinate of a particle).</li>
* <li>v: (per-DOF) This is the current velocity associated with the degree of
* freedom (the x, y, or z component of a particle's velocity).</li>
* <li>f: (per-DOF, read-only) This is the current force acting on the degree of
* freedom (the x, y, or z component of the force on a particle).</li>
* <li>m: (per-DOF, read-only) This is the mass of the particle the degree of
* freedom is associated with.</li>
* <li>uniform: (either global or per-DOF, read-only) This is a uniformly
* distributed random number between 0 and 1. Every time an expression is
* evaluated, a different value will be used. When used in a per-DOF
* expression, a different value will be used for every degree of freedom.
* Note, however, that if this variable appears multiple times in a single
* expression, the <i>same</i> value is used everywhere it appears in that
* expression.</li>
* <li>gaussian: (either global or per-DOF, read-only) This is a Gaussian
* distributed random number with mean 0 and variance 1. Every time an expression
* is evaluated, a different value will be used. When used in a per-DOF
* expression, a different value will be used for every degree of freedom.
* Note, however, that if this variable appears multiple times in a single
* expression, the <i>same</i> value is used everywhere it appears in that
* expression.</li>
* <li>A global variable is created for every adjustable parameter defined
* in the integrator's Context.</li>
* </ul>
*
* The following example uses a CustomIntegrator to implement a velocity Verlet
* integrator:
*
* <tt><pre>
* CustomIntegrator integrator;
* integrator.addComputePerDof("v", "v+0.5*dt*f/m");
* integrator.addComputePerDof("x", "x+dt*v");
* integrator.addComputePerDof("v", "v+0.5*dt*f/m");
* </pre></tt>
*
* The first step updates the velocities based on the current forces.
* The second step updates the positions based on the new velocities, and the
* third step updates the velocities again. Although the first and third steps
* look identical, the forces used in them are different. You do not need to
* tell the integrator that; it will recognize that the positions have changed
* and know to recompute the forces automatically.
*
* The above example has two problems. First, it does not respect distance
* constraints. To make the integrator work with constraints, you need to add
* extra steps to tell it when and how to apply them. Second, it never gives
* Forces an opportunity to update the context state. This should be done every
* time step so that, for example, an AndersenThermostat can randomize velocities
* or a MonteCarloBarostat can scale particle positions. You need to add a
* step to tell the integrator when to do this. The following example corrects
* both these problems:
*
* <tt><pre>
* CustomIntegrator integrator;
* integrator.addComputePerDof("v", "v+0.5*dt*f/m");
* integrator.addComputePerDof("x", "x+dt*v");
* integrator.addConstrainPositions();
* integrator.addUpdateContextState();
* integrator.addComputePerDof("v", "v+0.5*dt*f/m");
* integrator.addConstrainVelocities();
* </pre></tt>
*
* This integrator includes two steps that require forces (the two velocity
* updates) and three steps that can potentially change particle positions and
* thus invalidate the forces (the position update, applying position constraints,
* and allowing Forces to update the context state). We put all three of these
* steps together to minimize the number of force computations needed. If we had
* put addUpdateContextState() at the beginning of the algorithm instead, that would
* risk invalidating the forces just before the first velocity update, thus
* requiring two force evaluations per time step instead of one.
*
* Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
* functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, step. All trigonometric functions
* are defined in radians, and log is the natural logarithm. step(x) = 0 if x is less than 0, 1 otherwise. An expression
* may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
*/
class OPENMM_EXPORT CustomIntegrator : public Integrator {
public:
/**
* This is an enumeration of the different types of computations that may appear in an integration algorithm.
*/
enum ComputationType {
/**
* Compute an expression and store it in a global variable.
*/
ComputeGlobal = 0,
/**
* Compute an expression for every degree of freedom and store it in a per-DOF variable.
*/
ComputePerDof = 1,
/**
* Compute an expression for every degree of freedom, sum the values, and store the result in a global variable.
*/
ComputeSum = 2,
/**
* Update particle positions so all constraints are satisfied.
*/
ConstrainPositions = 3,
/**
* Update particle velocities so the net velocity along all constraints is 0.
*/
ConstrainVelocities = 4,
/**
* Allow Forces to update the context state.
*/
UpdateContextState = 5
};
/**
* Create a CustomIntegrator.
*
* @param stepSize the step size with which to integrate the system (in picoseconds)
*/
CustomIntegrator(double stepSize);
/**
* Get the number of global variables that have been defined.
*/
int getNumGlobalVariables() const {
return globalNames.size();
}
/**
* Get the number of per-DOF variables that have been defined.
*/
int getNumPerDofVariables() const {
return perDofNames.size();
}
/**
* Get the number of computation steps that have been added.
*/
int getNumComputations() const {
return computations.size();
}
/**
* Define a new global variable.
*
* @param name the name of the variable
* @param initialValue the variable will initially be set to this value
* @return the index of the variable that was added
*/
int addGlobalVariable(const std::string& name, double initialValue);
/**
* Get the name of a global variable.
*
* @param index the index of the variable to get
* @return the name of the variable
*/
std::string getGlobalVariableName(int index) const;
/**
* Define a new per-DOF variable.
*
* @param name the name of the variable
* @param initialValue the variable will initially be set to this value for
* all degrees of freedom
* @return the index of the variable that was added
*/
int addPerDofVariable(const std::string& name, double initialValue);
/**
* Get the name of a per-DOF variable.
*
* @param index the index of the variable to get
* @return the name of the variable
*/
std::string getPerDofVariableName(int index) const;
/**
* Get the current value of a global variable.
*
* @param index the index of the variable to get
* @return the current value of the variable
*/
double getGlobalVariable(int index) const;
/**
* Set the value of a global variable.
*
* @param index the index of the variable to get
* @param value the new value of the variable
*/
void setGlobalVariable(int index, double value);
/**
* Get the value of a per-DOF variable.
*
* @param index the index of the variable to get
* @param values the values of the variable for all degrees of freedom
* are stored into this
*/
void getPerDofVariable(int index, std::vector<Vec3>& values) const;
/**
* Set the value of a per-DOF variable.
*
* @param index the index of the variable to get
* @param values the new values of the variable for all degrees of freedom
*/
void setPerDofVariable(int index, const std::vector<Vec3>& values);
/**
* Add a step to the integration algorithm that computes a global value.
*
* @param variable the global variable to store the computed value into
* @param expression a mathematical expression involving only global variables.
* In each integration step, its value is computed and
* stored into the specified variable.
* @return the index of the step that was added
*/
int addComputeGlobal(const std::string& variable, const std::string& expression);
/**
* Add a step to the integration algorithm that computes a per-DOF value.
*
* @param variable the per-DOF variable to store the computed value into
* @param expression a mathematical expression involving both global and
* per-DOF variables. In each integration step, its value
* is computed for every degree of freedom and stored into
* the specified variable.
* @return the index of the step that was added
*/
int addComputePerDof(const std::string& variable, const std::string& expression);
/**
* Add a step to the integration algorithm that computes a sum over degrees of freedom.
*
* @param variable the global variable to store the computed value into
* @param expression a mathematical expression involving both global and
* per-DOF variables. In each integration step, its value
* is computed for every degree of freedom. Those values
* are then added together, and the sum is stored in the
* specified variable.
* @return the index of the step that was added
*/
int addComputeSum(const std::string& variable, const std::string& expression);
/**
* Add a step to the integration algorithm that updates particle positions so
* all constraints are satisfied.
*
* @return the index of the step that was added
*/
int addConstrainPositions();
/**
* Add a step to the integration algorithm that updates particle velocities
* so the net velocity along all constraints is 0.
*
* @return the index of the step that was added
*/
int addConstrainVelocities();
/**
* Add a step to the integration algorithm that allows Forces to update the
* context state.
*
* @return the index of the step that was added
*/
int addUpdateContextState();
/**
* Get the details of a computation step that has been added to the integration algorithm.
*
* @param index the index of the computation step to get
* @param type on exit, the type of computation this step performs
* @param variable on exit, the variable into which this step stores its result. If this
* step does not store a result in a variable, this will be an
* empty string.
* @param expression on exit, the expression this step evaluates. If this step does not
* evaluate an expression, this will be an empty string.
*/
void getComputationStep(int index, ComputationType& type, std::string& variable, std::string& expression) const;
/**
* Get the random number seed. See setRandomNumberSeed() for details.
*/
int getRandomNumberSeed() const {
return randomNumberSeed;
}
/**
* Set the random number seed. The precise meaning of this parameter is undefined, and is left up
* to each Platform to interpret in an appropriate way. It is guaranteed that if two simulations
* are run with different random number seeds, the sequence of random numbers 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;
}
/**
* Advance a simulation through time by taking a series of time steps.
*
* @param steps the number of time steps to take
*/
void step(int steps);
protected:
/**
* This will be called by the Context when it is created. It informs the Integrator
* of what context it will be integrating, and gives it a chance to do any necessary initialization.
* It will also get called again if the application calls reinitialize() on the Context.
*/
void initialize(ContextImpl& context);
/**
* When the user modifies the state, we need to mark that the forces need to be recalculated.
*/
void stateChanged(State::DataType changed);
/**
* Get the names of all Kernels used by this Integrator.
*/
std::vector<std::string> getKernelNames();
private:
class ComputationInfo;
std::vector<std::string> globalNames;
std::vector<std::string> perDofNames;
mutable std::vector<double> globalValues;
std::vector<std::vector<Vec3> > perDofValues;
std::vector<ComputationInfo> computations;
mutable bool globalsAreCurrent;
int randomNumberSeed;
bool forcesAreValid;
ContextImpl* context;
Context* owner;
Kernel kernel;
};
/**
* This is an internal class used to record information about a computation step.
* @private
*/
class CustomIntegrator::ComputationInfo {
public:
ComputationType type;
std::string variable, expression;
ComputationInfo() {
}
ComputationInfo(ComputationType type, const std::string& variable, const std::string& expression) :
type(type), variable(variable), expression(expression) {
}
};
} // namespace OpenMM
#endif /*OPENMM_CUSTOMINTEGRATOR_H_*/
openmmapi/include/openmm/internal/ContextImpl.h
View file @ dfba24ea
......@@ -118,6 +118,10 @@ public:
* @param forces on exit, this contains the forces
*/
void getForces(std::vector<Vec3>& forces);
/**
* Get the set of all adjustable parameters and their values
*/
const std::map<std::string, double>& getParameters() const;
/**
* Get the value of an adjustable parameter. If there is no parameter with the specified name, this
* throws an exception.
......
openmmapi/src/ContextImpl.cpp
View file @ dfba24ea
......@@ -123,6 +123,10 @@ void ContextImpl::getForces(std::vector<Vec3>& forces) {
dynamic_cast<UpdateStateDataKernel&>(updateStateDataKernel.getImpl()).getForces(*this, forces);
}
const std::map<std::string, double>& ContextImpl::getParameters() const {
return parameters;
}
double ContextImpl::getParameter(std::string name) {
if (parameters.find(name) == parameters.end())
throw OpenMMException("Called getParameter() with invalid parameter name");
......
openmmapi/src/CustomIntegrator.cpp 0 → 100644
View file @ dfba24ea
/* -------------------------------------------------------------------------- *
* 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) 2011 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/CustomIntegrator.h"
#include "openmm/Context.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/kernels.h"
#include <ctime>
#include <string>
using namespace OpenMM;
using std::string;
using std::vector;
CustomIntegrator::CustomIntegrator(double stepSize) : owner(NULL), globalsAreCurrent(true), forcesAreValid(false) {
setStepSize(stepSize);
setConstraintTolerance(1e-4);
setRandomNumberSeed((int) time(NULL));
}
void CustomIntegrator::initialize(ContextImpl& contextRef) {
if (owner != NULL && &contextRef.getOwner() != owner)
throw OpenMMException("This Integrator is already bound to a context");
context = &contextRef;
owner = &contextRef.getOwner();
kernel = context->getPlatform().createKernel(IntegrateCustomStepKernel::Name(), contextRef);
dynamic_cast<IntegrateCustomStepKernel&>(kernel.getImpl()).initialize(contextRef.getSystem(), *this);
dynamic_cast<IntegrateCustomStepKernel&>(kernel.getImpl()).setGlobalVariables(contextRef, globalValues);
for (int i = 0; i < (int) perDofValues.size(); i++)
dynamic_cast<IntegrateCustomStepKernel&>(kernel.getImpl()).setPerDofVariable(contextRef, i, perDofValues[i]);
}
void CustomIntegrator::stateChanged(State::DataType changed) {
forcesAreValid = false;
}
vector<string> CustomIntegrator::getKernelNames() {
std::vector<std::string> names;
names.push_back(IntegrateCustomStepKernel::Name());
return names;
}
void CustomIntegrator::step(int steps) {
globalsAreCurrent = false;
for (int i = 0; i < steps; ++i) {
context->updateContextState();
dynamic_cast<IntegrateCustomStepKernel&>(kernel.getImpl()).execute(*context, *this, forcesAreValid);
}
}
int CustomIntegrator::addGlobalVariable(const std::string& name, double initialValue) {
if (owner != NULL)
throw OpenMMException("The integrator cannot be modified after it is bound to a context");
globalNames.push_back(name);
globalValues.push_back(initialValue);
return globalNames.size()-1;
}
std::string CustomIntegrator::getGlobalVariableName(int index) const {
if (index < 0 || index >= globalNames.size())
throw OpenMMException("Index out of range");
return globalNames[index];
}
int CustomIntegrator::addPerDofVariable(const std::string& name, double initialValue) {
if (owner != NULL)
throw OpenMMException("The integrator cannot be modified after it is bound to a context");
perDofNames.push_back(name);
perDofValues.push_back(vector<Vec3>(1, Vec3(initialValue, initialValue, initialValue)));
return perDofNames.size()-1;
}
std::string CustomIntegrator::getPerDofVariableName(int index) const {
if (index < 0 || index >= perDofNames.size())
throw OpenMMException("Index out of range");
return perDofNames[index];
}
double CustomIntegrator::getGlobalVariable(int index) const {
if (index < 0 || index >= globalNames.size())
throw OpenMMException("Index out of range");
if (owner != NULL && !globalsAreCurrent) {
dynamic_cast<const IntegrateCustomStepKernel&>(kernel.getImpl()).getGlobalVariables(*context, globalValues);
globalsAreCurrent = true;
}
return globalValues[index];
}
void CustomIntegrator::setGlobalVariable(int index, double value) {
if (index < 0 || index >= globalNames.size())
throw OpenMMException("Index out of range");
if (owner != NULL && !globalsAreCurrent) {
dynamic_cast<IntegrateCustomStepKernel&>(kernel.getImpl()).getGlobalVariables(*context, globalValues);
globalsAreCurrent = true;
}
globalValues[index] = value;
dynamic_cast<IntegrateCustomStepKernel&>(kernel.getImpl()).setGlobalVariables(*context, globalValues);
}
void CustomIntegrator::getPerDofVariable(int index, std::vector<Vec3>& values) const {
if (index < 0 || index >= perDofNames.size())
throw OpenMMException("Index out of range");
if (owner == NULL)
values = perDofValues[index];
else
dynamic_cast<const IntegrateCustomStepKernel&>(kernel.getImpl()).getPerDofVariable(*context, index, values);
}
void CustomIntegrator::setPerDofVariable(int index, const std::vector<Vec3>& values) {
if (index < 0 || index >= perDofNames.size())
throw OpenMMException("Index out of range");
if (owner == NULL)
perDofValues[index] = values;
else
dynamic_cast<IntegrateCustomStepKernel&>(kernel.getImpl()).setPerDofVariable(*context, index, values);
}
int CustomIntegrator::addComputeGlobal(const std::string& variable, const std::string& expression) {
if (owner != NULL)
throw OpenMMException("The integrator cannot be modified after it is bound to a context");
computations.push_back(ComputationInfo(ComputeGlobal, variable, expression));
return computations.size()-1;
}
int CustomIntegrator::addComputePerDof(const std::string& variable, const std::string& expression) {
if (owner != NULL)
throw OpenMMException("The integrator cannot be modified after it is bound to a context");
computations.push_back(ComputationInfo(ComputePerDof, variable, expression));
return computations.size()-1;
}
int CustomIntegrator::addComputeSum(const std::string& variable, const std::string& expression) {
if (owner != NULL)
throw OpenMMException("The integrator cannot be modified after it is bound to a context");
computations.push_back(ComputationInfo(ComputeSum, variable, expression));
return computations.size()-1;
}
int CustomIntegrator::addConstrainPositions() {
if (owner != NULL)
throw OpenMMException("The integrator cannot be modified after it is bound to a context");
computations.push_back(ComputationInfo(ConstrainPositions, "", ""));
return computations.size()-1;
}
int CustomIntegrator::addConstrainVelocities() {
if (owner != NULL)
throw OpenMMException("The integrator cannot be modified after it is bound to a context");
computations.push_back(ComputationInfo(ConstrainVelocities, "", ""));
return computations.size()-1;
}
int CustomIntegrator::addUpdateContextState() {
if (owner != NULL)
throw OpenMMException("The integrator cannot be modified after it is bound to a context");
computations.push_back(ComputationInfo(UpdateContextState, "", ""));
return computations.size()-1;
}
void CustomIntegrator::getComputationStep(int index, ComputationType& type, std::string& variable, std::string& expression) const {
if (index < 0 || index >= computations.size())
throw OpenMMException("Index out of range");
type = computations[index].type;
variable = computations[index].variable;
expression = computations[index].expression;
}
platforms/reference/src/ReferenceKernelFactory.cpp
View file @ dfba24ea
......@@ -86,6 +86,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
return new ReferenceIntegrateVariableLangevinStepKernel(name, platform, data);
if (name == IntegrateVariableVerletStepKernel::Name())
return new ReferenceIntegrateVariableVerletStepKernel(name, platform, data);
if (name == IntegrateCustomStepKernel::Name())
return new ReferenceIntegrateCustomStepKernel(name, platform, data);
if (name == ApplyAndersenThermostatKernel::Name())
return new ReferenceApplyAndersenThermostatKernel(name, platform);
if (name == ApplyMonteCarloBarostatKernel::Name())
......
platforms/reference/src/ReferenceKernels.cpp
View file @ dfba24ea
......@@ -40,6 +40,7 @@
#include "SimTKReference/ReferenceCMAPTorsionIxn.h"
#include "SimTKReference/ReferenceCustomAngleIxn.h"
#include "SimTKReference/ReferenceCustomBondIxn.h"
#include "SimTKReference/ReferenceCustomDynamics.h"
#include "SimTKReference/ReferenceCustomExternalIxn.h"
#include "SimTKReference/ReferenceCustomGBIxn.h"
#include "SimTKReference/ReferenceCustomHbondIxn.h"
......@@ -1255,6 +1256,7 @@ void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const Ver
dynamics->setReferenceConstraintAlgorithm(constraints);
prevStepSize = stepSize;
}
constraints->setTolerance(integrator.getConstraintTolerance());
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
data.time += stepSize;
data.stepCount++;
......@@ -1318,6 +1320,7 @@ void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const L
prevFriction = friction;
prevStepSize = stepSize;
}
constraints->setTolerance(integrator.getConstraintTolerance());
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
data.time += stepSize;
data.stepCount++;
......@@ -1380,6 +1383,7 @@ void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const B
prevFriction = friction;
prevStepSize = stepSize;
}
constraints->setTolerance(integrator.getConstraintTolerance());
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
data.time += stepSize;
data.stepCount++;
......@@ -1439,6 +1443,7 @@ void ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context,
prevFriction = friction;
prevErrorTol = errorTol;
}
constraints->setTolerance(integrator.getConstraintTolerance());
RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
data.time += dynamics->getDeltaT();
......@@ -1495,6 +1500,7 @@ void ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, c
dynamics->setReferenceConstraintAlgorithm(constraints);
prevErrorTol = errorTol;
}
constraints->setTolerance(integrator.getConstraintTolerance());
RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
data.time += dynamics->getDeltaT();
......@@ -1503,6 +1509,93 @@ void ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, c
data.stepCount++;
}
ReferenceIntegrateCustomStepKernel::~ReferenceIntegrateCustomStepKernel() {
if (dynamics)
delete dynamics;
if (constraints)
delete constraints;
if (constraintIndices)
disposeIntArray(constraintIndices, numConstraints);
if (constraintDistances)
delete[] constraintDistances;
}
void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const CustomIntegrator& integrator) {
int numParticles = system.getNumParticles();
masses.resize(numParticles);
for (int i = 0; i < numParticles; ++i)
masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
numConstraints = system.getNumConstraints();
constraintIndices = allocateIntArray(numConstraints, 2);
constraintDistances = new RealOpenMM[numConstraints];
for (int i = 0; i < numConstraints; ++i) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
constraintIndices[i][0] = particle1;
constraintIndices[i][1] = particle2;
constraintDistances[i] = static_cast<RealOpenMM>(distance);
}
perDofValues.resize(integrator.getNumPerDofVariables());
for (int i = 0; i < (int) perDofValues.size(); i++)
perDofValues[i].resize(numParticles);
}
void ReferenceIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& velData = extractVelocities(context);
vector<RealVec>& forceData = extractForces(context);
if (dynamics == 0) {
// Create the computation objects.
dynamics = new ReferenceCustomDynamics(context.getSystem().getNumParticles(), integrator);
vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
findAnglesForCCMA(context.getSystem(), angles);
constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
dynamics->setReferenceConstraintAlgorithm(constraints);
}
// Record global variables.
map<string, double> globals;
globals["dt"] = integrator.getStepSize();
for (int i = 0; i < integrator.getNumGlobalVariables(); i++)
globals[integrator.getGlobalVariableName(i)] = globalValues[i];
// Execute the step.
constraints->setTolerance(integrator.getConstraintTolerance());
dynamics->update(context, context.getSystem().getNumParticles(), posData, velData, forceData, masses, globals, perDofValues, forcesAreValid);
// Record changed global variables.
integrator.setStepSize(globals["dt"]);
for (int i = 0; i < (int) globalValues.size(); i++)
globalValues[i] = globals[integrator.getGlobalVariableName(i)];
data.time += dynamics->getDeltaT();
data.stepCount++;
}
void ReferenceIntegrateCustomStepKernel::getGlobalVariables(ContextImpl& context, vector<double>& values) const {
values = globalValues;
}
void ReferenceIntegrateCustomStepKernel::setGlobalVariables(ContextImpl& context, const vector<double>& values) {
globalValues = values;
}
void ReferenceIntegrateCustomStepKernel::getPerDofVariable(ContextImpl& context, int variable, vector<Vec3>& values) const {
values.resize(perDofValues[variable].size());
for (int i = 0; i < (int) values.size(); i++)
values[i] = perDofValues[variable][i];
}
void ReferenceIntegrateCustomStepKernel::setPerDofVariable(ContextImpl& context, int variable, const vector<Vec3>& values) {
perDofValues[variable].resize(values.size());
for (int i = 0; i < (int) values.size(); i++)
perDofValues[variable][i] = values[i];
}
ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
if (thermostat)
delete thermostat;
......
platforms/reference/src/ReferenceKernels.h
View file @ dfba24ea
......@@ -49,6 +49,7 @@ class ReferenceMonteCarloBarostat;
class ReferenceVariableStochasticDynamics;
class ReferenceVariableVerletDynamics;
class ReferenceVerletDynamics;
class ReferenceCustomDynamics;
namespace OpenMM {
......@@ -857,6 +858,74 @@ private:
double prevErrorTol;
};
/**
* This kernel is invoked by CustomIntegrator to take one time step.
*/
class ReferenceIntegrateCustomStepKernel : public IntegrateCustomStepKernel {
public:
ReferenceIntegrateCustomStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateCustomStepKernel(name, platform),
data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
}
~ReferenceIntegrateCustomStepKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param integrator the CustomIntegrator this kernel will be used for
*/
void initialize(const System& system, const CustomIntegrator& integrator);
/**
* Execute the kernel.
*
* @param context the context in which to execute this kernel
* @param integrator the CustomIntegrator this kernel is being used for
* @param forcesAreValid if the context has been modified since the last time step, this will be
* false to show that cached forces are invalid and must be recalculated.
* On exit, this should specify whether the cached forces are valid at the
* end of the step.
*/
void execute(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid);
/**
* Get the values of all global variables.
*
* @param context the context in which to execute this kernel
* @param values on exit, this contains the values
*/
void getGlobalVariables(ContextImpl& context, std::vector<double>& values) const;
/**
* Set the values of all global variables.
*
* @param context the context in which to execute this kernel
* @param values a vector containing the values
*/
void setGlobalVariables(ContextImpl& context, const std::vector<double>& values);
/**
* Get the values of a per-DOF variable.
*
* @param context the context in which to execute this kernel
* @param variable the index of the variable to get
* @param values on exit, this contains the values
*/
void getPerDofVariable(ContextImpl& context, int variable, std::vector<Vec3>& values) const;
/**
* Set the values of a per-DOF variable.
*
* @param context the context in which to execute this kernel
* @param variable the index of the variable to get
* @param values a vector containing the values
*/
void setPerDofVariable(ContextImpl& context, int variable, const std::vector<Vec3>& values);
private:
ReferencePlatform::PlatformData& data;
ReferenceCustomDynamics* dynamics;
ReferenceConstraintAlgorithm* constraints;
std::vector<RealOpenMM> masses, globalValues;
std::vector<std::vector<OpenMM::RealVec> > perDofValues;
RealOpenMM* constraintDistances;
int** constraintIndices;
int numConstraints;
};
/**
* This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
*/
......
platforms/reference/src/ReferencePlatform.cpp
View file @ dfba24ea
......@@ -65,6 +65,7 @@ ReferencePlatform::ReferencePlatform() {
registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);
registerKernelFactory(IntegrateVariableVerletStepKernel::Name(), factory);
registerKernelFactory(IntegrateCustomStepKernel::Name(), factory);
registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
registerKernelFactory(ApplyMonteCarloBarostatKernel::Name(), factory);
registerKernelFactory(CalcKineticEnergyKernel::Name(), factory);
......
platforms/reference/src/SimTKReference/ReferenceCustomDynamics.cpp 0 → 100644
View file @ dfba24ea
/* Portions copyright (c) 2011 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* 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 "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "openmm/OpenMMException.h"
#include "lepton/Parser.h"
#include "ReferenceCustomDynamics.h"
#include "lepton/ParsedExpression.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/ForceImpl.h"
#include "lepton/Operation.h"
#include <set>
using namespace std;
using namespace OpenMM;
/**---------------------------------------------------------------------------------------
ReferenceCustomDynamics constructor
@param numberOfAtoms number of atoms
@param integrator the integrator definition to use
--------------------------------------------------------------------------------------- */
ReferenceCustomDynamics::ReferenceCustomDynamics(int numberOfAtoms, const CustomIntegrator& integrator) :
ReferenceDynamics(numberOfAtoms, integrator.getStepSize(), 0.0), integrator(integrator) {
sumBuffer.resize(numberOfAtoms);
stepType.resize(integrator.getNumComputations());
stepVariable.resize(integrator.getNumComputations());
stepExpression.resize(integrator.getNumComputations());
for (int i = 0; i < integrator.getNumComputations(); i++) {
string expression;
integrator.getComputationStep(i, stepType[i], stepVariable[i], expression);
if (expression.length() > 0)
stepExpression[i] = Lepton::Parser::parse(expression).createProgram();
}
}
/**---------------------------------------------------------------------------------------
ReferenceCustomDynamics destructor
--------------------------------------------------------------------------------------- */
ReferenceCustomDynamics::~ReferenceCustomDynamics() {
}
/**---------------------------------------------------------------------------------------
Update -- driver routine for performing Custom dynamics update of coordinates
and velocities
@param context the context this integrator is updating
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param globals a map containing values of global variables
@param forcesAreValid whether the current forces are valid or need to be recomputed
--------------------------------------------------------------------------------------- */
void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses,
map<string, RealOpenMM>& globals, vector<vector<RealVec> >& perDof, bool& forcesAreValid){
int numSteps = stepType.size();
globals.insert(context.getParameters().begin(), context.getParameters().end());
if (invalidatesForces.size() == 0) {
// The first time this is called, work out when to recompute forces and energy. First build a
// list of every step that invalidates the forces.
invalidatesForces.resize(numSteps, false);
needsForces.resize(numSteps, false);
needsEnergy.resize(numSteps, false);
set<string> affectsForce;
affectsForce.insert("x");
for (vector<ForceImpl*>::const_iterator iter = context.getForceImpls().begin(); iter != context.getForceImpls().end(); ++iter) {
const map<string, double> params = (*iter)->getDefaultParameters();
for (map<string, double>::const_iterator param = params.begin(); param != params.end(); ++param)
affectsForce.insert(param->first);
}
for (int i = 0; i < numSteps; i++)
invalidatesForces[i] = (stepType[i] == CustomIntegrator::ConstrainPositions || affectsForce.find(stepVariable[i]) != affectsForce.end());
// Make a list of which steps require valid forces or energy to be known.
for (int i = 0; i < numSteps; i++) {
if (stepType[i] == CustomIntegrator::ComputeGlobal || stepType[i] == CustomIntegrator::ComputePerDof || stepType[i] == CustomIntegrator::ComputeSum) {
for (int j = 0; j < stepExpression[i].getNumOperations(); j++) {
const Lepton::Operation& op = stepExpression[i].getOperation(j);
if (op.getId() == Lepton::Operation::VARIABLE) {
if (op.getName() == "f")
needsForces[i] = true;
else if (op.getName() == "energy")
needsEnergy[i] = true;
}
}
}
}
// Build the list of inverse masses.
inverseMasses.resize(numberOfAtoms);
for (int i = 0; i < numberOfAtoms; i++)
inverseMasses[i] = 1.0/masses[i];
}
// Loop over steps and execute them.
for (int i = 0; i < numSteps; i++) {
if ((needsForces[i] || needsEnergy[i]) && !forcesAreValid) {
// Recompute forces and or energy. Figure out what is actually needed
// between now and the next time they get invalidated again.
bool computeForce = false, computeEnergy = false;
for (int j = i; ; j++) {
if (needsForces[j])
computeForce = true;
if (needsEnergy[j])
computeEnergy = true;
if (invalidatesForces[j])
break;
if (j == numSteps-1)
j = -1;
if (j == i-1)
break;
}
recordChangedParameters(context, globals);
RealOpenMM e = context.calcForcesAndEnergy(computeForce, computeEnergy);
if (computeEnergy)
energy = e;
forcesAreValid = true;
}
globals["energy"] = energy;
// Execute the step.
switch (stepType[i]) {
case CustomIntegrator::ComputeGlobal: {
map<string, RealOpenMM> variables = globals;
variables["uniform"] = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber();
variables["gaussian"] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
globals[stepVariable[i]] = stepExpression[i].evaluate(variables);
break;
}
case CustomIntegrator::ComputePerDof: {
vector<RealVec>* results = NULL;
if (stepVariable[i] == "x")
results = &atomCoordinates;
else if (stepVariable[i] == "v")
results = &velocities;
else {
for (int j = 0; j < integrator.getNumPerDofVariables(); j++)
if (stepVariable[i] == integrator.getPerDofVariableName(j))
results = &perDof[j];
}
if (results == NULL)
throw OpenMMException("Illegal per-DOF output variable: "+stepVariable[i]);
computePerDof(numberOfAtoms, *results, atomCoordinates, velocities, forces, masses, globals, perDof, stepExpression[i]);
break;
}
case CustomIntegrator::ComputeSum: {
computePerDof(numberOfAtoms, sumBuffer, atomCoordinates, velocities, forces, masses, globals, perDof, stepExpression[i]);
RealOpenMM sum = 0.0;
for (int j = 0; j < numberOfAtoms; j++)
sum += sumBuffer[j][0]+sumBuffer[j][1]+sumBuffer[j][2];
globals[stepVariable[i]] = sum;
break;
}
case CustomIntegrator::ConstrainPositions: {
getReferenceConstraintAlgorithm()->apply(numberOfAtoms, atomCoordinates, atomCoordinates, inverseMasses);
}
case CustomIntegrator::UpdateContextState: {
recordChangedParameters(context, globals);
context.updateContextState();
globals.insert(context.getParameters().begin(), context.getParameters().end());
}
}
if (invalidatesForces[i])
forcesAreValid = false;
}
incrementTimeStep();
recordChangedParameters(context, globals);
}
RealOpenMM ReferenceCustomDynamics::computePerDof(int numberOfAtoms, vector<RealVec>& results, const vector<RealVec>& atomCoordinates,
const vector<RealVec>& velocities, const vector<RealVec>& forces, const vector<RealOpenMM>& masses,
const map<string, RealOpenMM>& globals, const vector<vector<RealVec> >& perDof, const Lepton::ExpressionProgram& expression) {
// Loop over all degrees of freedom.
map<string, RealOpenMM> variables = globals;
for (int i = 0; i < numberOfAtoms; i++) {
variables["m"] = masses[i];
for (int j = 0; j < 3; j++) {
// Compute the expression.
variables["x"] = atomCoordinates[i][j];
variables["v"] = velocities[i][j];
variables["f"] = forces[i][j];
variables["uniform"] = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber();
variables["gaussian"] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
for (int k = 0; k < (int) perDof.size(); k++)
variables[integrator.getPerDofVariableName(k)] = perDof[k][i][j];
results[i][j] = expression.evaluate(variables);
}
}
}
/**
* Check which context parameters have changed and register them with the context.
*/
void ReferenceCustomDynamics::recordChangedParameters(OpenMM::ContextImpl& context, std::map<std::string, RealOpenMM>& globals) {
for (map<string, double>::const_iterator iter = context.getParameters().begin(); iter != context.getParameters().end(); ++iter) {
string name = iter->first;
double value = globals[name];
if (value != iter->second)
context.setParameter(name, globals[name]);
}
}
platforms/reference/src/SimTKReference/ReferenceCustomDynamics.h 0 → 100644
View file @ dfba24ea
/* Portions copyright (c) 2011 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* 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.
*/
#ifndef __ReferenceCustomDynamics_H__
#define __ReferenceCustomDynamics_H__
#include "ReferenceDynamics.h"
#include "openmm/CustomIntegrator.h"
#include "openmm/internal/ContextImpl.h"
#include "lepton/ExpressionProgram.h"
#include <map>
#include <string>
#include <vector>
// ---------------------------------------------------------------------------------------
class ReferenceCustomDynamics : public ReferenceDynamics {
private:
const OpenMM::CustomIntegrator& integrator;
std::vector<RealOpenMM> inverseMasses;
std::vector<OpenMM::RealVec> sumBuffer;
std::vector<OpenMM::CustomIntegrator::ComputationType> stepType;
std::vector<std::string> stepVariable;
std::vector<Lepton::ExpressionProgram> stepExpression;
std::vector<bool> invalidatesForces, needsForces, needsEnergy;
RealOpenMM energy;
RealOpenMM computePerDof(int numberOfAtoms, std::vector<OpenMM::RealVec>& results, const std::vector<OpenMM::RealVec>& atomCoordinates,
const std::vector<OpenMM::RealVec>& velocities, const std::vector<OpenMM::RealVec>& forces, const std::vector<RealOpenMM>& masses,
const std::map<std::string, RealOpenMM>& globals, const std::vector<std::vector<OpenMM::RealVec> >& perDof, const Lepton::ExpressionProgram& expression);
void recordChangedParameters(OpenMM::ContextImpl& context, std::map<std::string, RealOpenMM>& globals);
public:
/**---------------------------------------------------------------------------------------
Constructor
@param numberOfAtoms number of atoms
@param integrator the integrator definition to use
--------------------------------------------------------------------------------------- */
ReferenceCustomDynamics(int numberOfAtoms, const OpenMM::CustomIntegrator& integrator);
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceCustomDynamics();
/**---------------------------------------------------------------------------------------
Update
@param context the context this integrator is updating
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param globals a map containing values of global variables
@param perDof the values of per-DOF variables
@param forcesAreValid whether the current forces are valid or need to be recomputed
--------------------------------------------------------------------------------------- */
void update(OpenMM::ContextImpl& context, int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses,
std::map<std::string, RealOpenMM>& globals, std::vector<std::vector<OpenMM::RealVec> >& perDof, bool& forcesAreValid);
};
// ---------------------------------------------------------------------------------------
#endif // __ReferenceCustomDynamics_H__
platforms/reference/tests/TestReferenceCustomIntegrator.cpp 0 → 100644
View file @ dfba24ea
/* -------------------------------------------------------------------------- *
* 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-2011 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 CustomIntegrator.
*/
#include "../../../tests/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/AndersenThermostat.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/CustomIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
/**
* Test a simple leapfrog integrator on a single bond.
*/
void testSingleBond() {
ReferencePlatform platform;
System system;
system.addParticle(2.0);
system.addParticle(2.0);
CustomIntegrator integrator(0.01);
integrator.addComputePerDof("v", "v+dt*f/m");
integrator.addComputePerDof("x", "x+dt*v");
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
// This is simply a harmonic oscillator, so compare it to the analytical solution.
const double freq = 1.0;;
State state = context.getState(State::Energy);
const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
for (int i = 0; i < 1000; ++i) {
state = context.getState(State::Positions | State::Velocities | State::Energy);
double time = state.getTime();
double expectedDist = 1.5+0.5*std::cos(freq*time);
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedDist, 0, 0), state.getPositions()[0], 0.02);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedDist, 0, 0), state.getPositions()[1], 0.02);
double expectedSpeed = -0.5*freq*std::sin(freq*time);
ASSERT_EQUAL_VEC(Vec3(-0.5*expectedSpeed, 0, 0), state.getVelocities()[0], 0.02);
ASSERT_EQUAL_VEC(Vec3(0.5*expectedSpeed, 0, 0), state.getVelocities()[1], 0.02);
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(1);
}
}
/**
* Test an integrator that enforces constraints.
*/
void testConstraints() {
const int numParticles = 8;
const double temp = 500.0;
ReferencePlatform platform;
System system;
CustomIntegrator integrator(0.002);
integrator.addPerDofVariable("oldx", 0);
integrator.addComputePerDof("v", "v+dt*f/m");
integrator.addComputePerDof("oldx", "x");
integrator.addComputePerDof("x", "x+dt*v");
integrator.addConstrainPositions();
integrator.addComputePerDof("v", "(x-oldx)/dt");
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
system.addParticle(i%2 == 0 ? 5.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
}
for (int i = 0; i < numParticles-1; ++i)
system.addConstraint(i, i+1, 1.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; ++i) {
positions[i] = Vec3(i/2, (i+1)/2, 0);
velocities[i] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Energy);
for (int j = 0; j < system.getNumConstraints(); ++j) {
int particle1, particle2;
double distance;
system.getConstraintParameters(j, particle1, particle2, distance);
Vec3 p1 = state.getPositions()[particle1];
Vec3 p2 = state.getPositions()[particle2];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
ASSERT_EQUAL_TOL(distance, dist, 2e-5);
}
double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1)
initialEnergy = energy;
else if (i > 1)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
integrator.step(1);
}
}
/**
* Test an integrator with an AndersenThermostat to see if updateContextState()
* is being handled correctly.
*/
void testWithThermostat() {
const int numParticles = 8;
const double temp = 100.0;
const double collisionFreq = 10.0;
const int numSteps = 10000;
ReferencePlatform platform;
System system;
CustomIntegrator integrator(0.005);
integrator.addUpdateContextState();
integrator.addComputePerDof("v", "v+dt*f/m");
integrator.addComputePerDof("x", "x+dt*v");
NonbondedForce* forceField = new NonbondedForce();
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);
AndersenThermostat* thermostat = new AndersenThermostat(temp, collisionFreq);
system.addForce(thermostat);
Context context(system, integrator, platform);
vector<Vec3> positions(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));
context.setPositions(positions);
// Let it equilibrate.
integrator.step(10000);
// Now run it for a while and see if the temperature is correct.
double ke = 0.0;
for (int i = 0; i < numSteps; ++i) {
State state = context.getState(State::Energy);
ke += state.getKineticEnergy();
integrator.step(1);
}
ke /= numSteps;
double expected = 0.5*numParticles*3*BOLTZ*temp;
ASSERT_USUALLY_EQUAL_TOL(expected, ke, 6/std::sqrt((double) numSteps));
}
/**
* Test a Monte Carlo integrator that uses global variables and depends on energy.
*/
void testMonteCarlo() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
CustomIntegrator integrator(0.1);
const double kT = BOLTZ*300.0;
integrator.addGlobalVariable("kT", kT);
integrator.addGlobalVariable("oldE", 0);
integrator.addGlobalVariable("accept", 0);
integrator.addPerDofVariable("oldx", 0);
integrator.addComputeGlobal("oldE", "energy");
integrator.addComputePerDof("oldx", "x");
integrator.addComputePerDof("x", "x+dt*gaussian");
integrator.addComputeGlobal("accept", "step(exp((oldE-energy)/kT)-uniform)");
integrator.addComputePerDof("x", "accept*x + (1-accept)*oldx");
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 2.0, 10.0);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
// Compute the histogram of distances and see if it satisfies a Boltzmann distribution.
const int numBins = 100;
const double maxDist = 4.0;
const int numIterations = 5000;
vector<int> counts(numBins, 0);
for (int i = 0; i < numIterations; ++i) {
integrator.step(10);
State state = context.getState(State::Positions);
Vec3 delta = state.getPositions()[0]-state.getPositions()[1];
double dist = sqrt(delta.dot(delta));
if (dist < maxDist)
counts[(int) (numBins*dist/maxDist)]++;
}
vector<double> expected(numBins, 0);
double sum = 0;
for (int i = 0; i < numBins; i++) {
double dist = (i+0.5)*maxDist/numBins;
expected[i] = dist*dist*exp(-5.0*(dist-2)*(dist-2)/kT);
sum += expected[i];
}
for (int i = 0; i < numBins; i++)
ASSERT_USUALLY_EQUAL_TOL((double) counts[i]/numIterations, expected[i]/sum, 0.01);
}
int main() {
try {
testSingleBond();
testConstraints();
testWithThermostat();
testMonteCarlo();
}
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