#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-2020 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 "TabulatedFunction.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: global variables
* have a single value, while per-DOF 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>
*
* Like all integrators, CustomIntegrator ignores any particle whose mass is 0.
* It is skipped when doing per-DOF computations, and is not included when
* computing sums over degrees of freedom.
*
* 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>energy0, energy1, energy2, ...: (global, read-only) This is similar to
* energy, but includes only the contribution from forces in one force group.
* A single computation step may only depend on a single energy variable
* (energy, energy0, energy1, etc.).</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>f0, f1, f2, ...: (per-DOF, read-only) This is similar to f, but includes
* only the contribution from forces in one force group. A single computation
* step may only depend on a single force variable (f, f0, f1, etc.).</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 same 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 same 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(0.001);
* 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, using the RATTLE algorithm to apply constraints:
*
* <tt><pre>
* CustomIntegrator integrator(0.001);
* integrator.addPerDofVariable("x1", 0);
* integrator.addUpdateContextState();
* integrator.addComputePerDof("v", "v+0.5*dt*f/m");
* integrator.addComputePerDof("x", "x+dt*v");
* integrator.addComputePerDof("x1", "x");
* integrator.addConstrainPositions();
* integrator.addComputePerDof("v", "v+0.5*dt*f/m+(x-x1)/dt");
* integrator.addConstrainVelocities();
* </pre></tt>
*
* CustomIntegrator can be used to implement multiple time step integrators. The
* following example shows an r-RESPA integrator. It assumes the quickly changing
* forces are in force group 0 and the slowly changing ones are in force group 1.
* It evaluates the "fast" forces four times as often as the "slow" forces.
*
* <tt><pre>
* CustomIntegrator integrator(0.004);
* integrator.addComputePerDof("v", "v+0.5*dt*f1/m");
* for (int i = 0; i < 4; i++) {
* integrator.addComputePerDof("v", "v+0.5*(dt/4)*f0/m");
* integrator.addComputePerDof("x", "x+(dt/4)*v");
* integrator.addComputePerDof("v", "v+0.5*(dt/4)*f0/m");
* }
* integrator.addComputePerDof("v", "v+0.5*dt*f1/m");
* </pre></tt>
*
* The sequence of computations in a CustomIntegrator can include flow control in
* the form of "if" and "while" blocks. The computations inside an "if" block
* are executed either zero or one times, depending on whether a condition is
* true. The computations inside a "while" block are executed repeatedly for as
* long as the condition remains true. Be very careful when writing "while"
* blocks; there is nothing to stop you from creating an infinite loop!
*
* For example, suppose you are writing a Monte Carlo algorithm. Assume you have
* already computed a new set of particle coordinates "xnew" and a step acceptance
* probability "acceptanceProbability". The following lines use an "if" block
* to decide whether to accept the step, and if it is accepted, store the new
* positions into "x".
*
* <tt><pre>
* integrator.beginIfBlock("uniform < acceptanceProbability");
* integrator.addComputePerDof("x", "xnew");
* integrator.endBlock();
* </pre></tt>
*
* The condition in an "if" or "while" block is evaluated globally, so it may
* only involve global variables, not per-DOF ones. It may use any of the
* following comparison operators: =, <. >, !=, <=, >=. Blocks may be nested
* inside each other.
*
* "Per-DOF" computations can also be thought of as per-particle computations
* that operate on three component vectors. For example, "x+dt*v" means to take
* the particle's velocity (a vector), multiply it by the step size, and add the
* position (also a vector). The result is a new vector that can be stored into
* a per-DOF variable with addComputePerDof(), or it can be summed over all
* components of all particles with addComputeSum(). Because the calculation is
* done on vectors, you can use functions that operate explicitly on vectors
* rather than just computing each component independently. For example, the
* following line uses a cross product to compute the angular momentum of each
* particle and stores it into a per-DOF variable.
*
* <tt><pre>
* integrator.addComputePerDof("angularMomentum", "m*cross(x, v)");
* </pre></tt>
*
* Here are two more examples that may be useful as starting points for writing
* your own integrators. The first one implements the algorithm used by the
* standard VerletIntegrator class. This is a leapfrog algorithm, in contrast
* to the velocity Verlet algorithm shown above, so it only requires applying
* constraints once in each time step.
*
* <tt><pre>
* CustomIntegrator integrator(dt);
* integrator.addPerDofVariable("x0", 0);
* integrator.addUpdateContextState();
* integrator.addComputePerDof("x0", "x");
* integrator.addComputePerDof("v", "v+dt*f/m");
* integrator.addComputePerDof("x", "x+dt*v");
* integrator.addConstrainPositions();
* integrator.addComputePerDof("v", "(x-x0)/dt");
* </pre></tt>
*
* The second one implements the algorithm used by the standard
* LangevinMiddleIntegrator class. kB is Boltzmann's constant.
*
* <tt><pre>
* CustomIntegrator integrator(dt);
* integrator.addGlobalVariable("a", exp(-friction*dt));
* integrator.addGlobalVariable("b", sqrt(1-exp(-2*friction*dt)));
* integrator.addGlobalVariable("kT", kB*temperature);
* integrator.addPerDofVariable("x1", 0);
* integrator.addUpdateContextState();
* integrator.addComputePerDof("v", "v + dt*f/m");
* integrator.addConstrainVelocities();
* integrator.addComputePerDof("x", "x + 0.5*dt*v");
* integrator.addComputePerDof("v", "a*v + b*sqrt(kT/m)*gaussian");
* integrator.addComputePerDof("x", "x + 0.5*dt*v");
* integrator.addComputePerDof("x1", "x");
* integrator.addConstrainPositions();
* integrator.addComputePerDof("v", "v + (x-x1)/dt");
* </pre></tt>
*
* Another feature of CustomIntegrator is that it can use derivatives of the
* potential energy with respect to context parameters. These derivatives are
* typically computed by custom forces, and are only computed if a Force object
* has been specifically told to compute them by calling addEnergyParameterDerivative()
* on it. CustomIntegrator provides a deriv() function for accessing these
* derivatives in global or per-DOF expressions. For example, "deriv(energy, lambda)"
* is the derivative of the total potentially energy with respect to the parameter
* lambda. You can also restrict it to a single force group by specifying a different
* variable for the first argument, such as "deriv(energy1, lambda)".
*
* An Integrator has one other job in addition to evolving the equations of motion:
* it defines how to compute the kinetic energy of the system. Depending on the
* integration method used, simply summing (mv^2)/2 over all degrees of
* freedom may not give the correct answer. For example, in a leapfrog integrator
* the velocities are "delayed" by half a time step, so the above formula would
* give the kinetic energy half a time step ago, not at the current time.
*
* Call setKineticEnergyExpression() to set an expression for the kinetic energy.
* It is computed for every degree of freedom (excluding ones whose mass is 0) and
* the result is summed. The default expression is "m*v*v/2", which is correct
* for many integrators.
*
* As example, the following line defines the correct way to compute kinetic energy
* when using a leapfrog algorithm:
*
* <tt><pre>
* integrator.setKineticEnergyExpression("m*v1*v1/2; v1=v+0.5*dt*f/m");
* </pre></tt>
*
* The kinetic energy expression may depend on the following pre-defined variables:
* x, v, f, m, dt. It also may depend on user-defined global and per-DOF variables,
* and on the values of adjustable parameters defined in the integrator's Context.
* It may not depend on any other variable, such as the potential energy,
* the force from a single force group, or a random number.
*
* 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, atan2, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select. All trigonometric functions
* are defined in radians, and log is the natural logarithm. step(x) = 0 if x is less than 0, 1 otherwise. delta(x) = 1 if x is 0, 0 otherwise.
* select(x,y,z) = z if x = 0, y otherwise. An expression may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
*
* Expressions used in ComputePerDof and ComputeSum steps can also use the following
* functions that operate on vectors: cross(a, b) is the cross product of two
* vectors; dot(a, b) is the dot product of two vectors; _x(a), _y(a), and _z(a)
* extract a single component from a vector; and vector(a, b, c) creates a new
* vector with the x component of the first argument, the y component of the
* second argument, and the z component of the third argument. Remember that every
* quantity appearing in a vector expression is a vector. Functions that appear
* to return a scalar really return a vector whose components are all the same.
* For example, _z(a) returns the vector (a.z, a.z, a.z). Likewise, wherever a
* constant appears in the expression, it really means a vector whose components
* all have the same value.
*
* In addition, you can call addTabulatedFunction() to define a new function based on tabulated values. You specify the function by
* creating a TabulatedFunction object. That function can then appear in expressions.
*/
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,
/**
* Begin an "if" block.
*/
IfBlockStart = 6,
/**
* Begin a while" block.
*/
WhileBlockStart = 7,
/**
* End an "if" or "while" block.
*/
BlockEnd = 8
};
/**
* Create a CustomIntegrator.
*
* @param stepSize the step size with which to integrate the system (in picoseconds)
*/
CustomIntegrator(double stepSize);
~CustomIntegrator();
/**
* 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();
}
/**
* Get the number of tabulated functions that have been defined.
*/
int getNumTabulatedFunctions() const {
return functions.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
*/
const 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
*/
const 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;
/**
* Get the current value of a global variable, specified by name.
*
* @param name the name of the variable to get
* @return the current value of the parameter
*/
double getGlobalVariableByName(const std::string& name) const;
/**
* Set the value of a global variable.
*
* @param index the index of the variable to set
* @param value the new value of the variable
*/
void setGlobalVariable(int index, double value);
/**
* Set the value of a global variable, specified by name.
*
* @param name the name of the variable to set
* @param value the new value of the variable
*/
void setGlobalVariableByName(const std::string& name, 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;
/**
* Get the value of a per-DOF variable, specified by name.
*
* @param name the name of the variable to get
* @param[out] values the values of the variable for all degrees of freedom
* are stored into this
*/
void getPerDofVariableByName(const std::string& name, std::vector<Vec3>& values) const;
/**
* Set the value of a per-DOF variable.
*
* @param index the index of the variable to set
* @param values the new values of the variable for all degrees of freedom
*/
void setPerDofVariable(int index, const std::vector<Vec3>& values);
/**
* Set the value of a per-DOF variable, specified by name.
*
* @param name the name of the variable to set
* @param values the new values of the variable for all degrees of freedom
*/
void setPerDofVariableByName(const std::string& name, 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();
/**
* Add a step which begins a new "if" block.
*
* @param condition a mathematical expression involving a comparison operator
* and global variables. All steps between this one and
* the end of the block are executed only if the condition
* is true.
*
* @return the index of the step that was added
*/
int beginIfBlock(const std::string& condition);
/**
* Add a step which begins a new "while" block.
*
* @param condition a mathematical expression involving a comparison operator
* and global variables. All steps between this one and
* the end of the block are executed repeatedly as long as
* the condition remains true.
*
* @return the index of the step that was added
*/
int beginWhileBlock(const std::string& condition);
/**
* Add a step which marks the end of the most recently begun "if" or "while"
* block.
*
* @return the index of the step that was added
*/
int endBlock();
/**
* 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[out] type the type of computation this step performs
* @param[out] variable 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[out] expression 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;
/**
* Add a tabulated function that may appear in expressions.
*
* @param name the name of the function as it appears in expressions
* @param function a TabulatedFunction object defining the function. The TabulatedFunction
* should have been created on the heap with the "new" operator. The
* integrator takes over ownership of it, and deletes it when the integrator itself is deleted.
* @return the index of the function that was added
*/
int addTabulatedFunction(const std::string& name, TabulatedFunction* function);
/**
* Get a const reference to a tabulated function that may appear in expressions.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
const TabulatedFunction& getTabulatedFunction(int index) const;
/**
* Get a reference to a tabulated function that may appear in expressions.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
TabulatedFunction& getTabulatedFunction(int index);
/**
* Get the name of a tabulated function that may appear in expressions.
*
* @param index the index of the function to get
* @return the name of the function as it appears in expressions
*/
const std::string& getTabulatedFunctionName(int index) const;
/**
* Get the expression to use for computing the kinetic energy. The expression is evaluated
* for every degree of freedom. Those values are then added together, and the sum
* is reported as the current kinetic energy.
*/
const std::string& getKineticEnergyExpression() const;
/**
* Set the expression to use for computing the kinetic energy. The expression is evaluated
* for every degree of freedom. Those values are then added together, and the sum
* is reported as the current kinetic energy.
*/
void setKineticEnergyExpression(const std::string& expression);
/**
* 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.
*
* 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;
}
/**
* 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);
/**
* This will be called by the Context when it is destroyed to let the Integrator do any necessary
* cleanup. It will also get called again if the application calls reinitialize() on the Context.
*/
void cleanup();
/**
* 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();
/**
* Compute the kinetic energy of the system at the current time.
*/
double computeKineticEnergy();
/**
* Get whether computeKineticEnergy() expects forces to have been computed.
*/
bool kineticEnergyRequiresForce() const;
/**
* This is called while writing checkpoints. It gives the integrator a chance to write
* its own data.
*/
void createCheckpoint(std::ostream& stream) const;
/**
* This is called while loading a checkpoint. The integrator should read in whatever
* data it wrote in createCheckpoint() and update its internal state accordingly.
*/
void loadCheckpoint(std::istream& stream);
/**
* This is called while creating a State. The Integrator should store the values
* of all time-varying parameters into the SerializationNode so they can be saved
* as part of the state.
*/
void serializeParameters(SerializationNode& node) const;
/**
* This is called when loading a previously saved State. The Integrator should
* load the values of all time-varying parameters from the SerializationNode. If
* the node contains parameters that are not defined for this Integrator, it should
* throw an exception.
*/
void deserializeParameters(const SerializationNode& node);
private:
class ComputationInfo;
class FunctionInfo;
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;
std::vector<FunctionInfo> functions;
std::string kineticEnergy;
mutable bool globalsAreCurrent;
int randomNumberSeed;
bool forcesAreValid, keNeedsForce;
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) {
}
};
/**
* This is an internal class used to record information about a tabulated function.
* @private
*/
class CustomIntegrator::FunctionInfo {
public:
std::string name;
TabulatedFunction* function;
FunctionInfo() {
}
FunctionInfo(const std::string& name, TabulatedFunction* function) : name(name), function(function) {
}
};
} // namespace OpenMM
#endif /*OPENMM_CUSTOMINTEGRATOR_H_*/