Commit 1b705ba9 authored by Peter Eastman's avatar Peter Eastman
Browse files

Minor cleanups to VariableVerletIntegrator

parent 09786e4b
...@@ -46,9 +46,9 @@ namespace OpenMM { ...@@ -46,9 +46,9 @@ namespace OpenMM {
* specified tolerance. This allows it to take larger steps on average than a fixed step size * specified tolerance. This allows it to take larger steps on average than a fixed step size
* integrator, while still maintaining comparable accuracy and stability. * integrator, while still maintaining comparable accuracy and stability.
* *
* It is best not to think of the value of the accuracy parameter as having any absolute meaning. * It is best not to think of the error tolerance as having any absolute meaning. It is just an
* It is just an adjustable parameter that affects the step size and integration accuracy. You * adjustable parameter that affects the step size and integration accuracy. You
* should try different values to find the largest value that produces a trajectory sufficiently * should try different values to find the largest onethat produces a trajectory sufficiently
* accurate for your purposes. 0.001 is often a good starting point. * accurate for your purposes. 0.001 is often a good starting point.
* *
* Unlike a fixed step size Verlet integrator, variable step size Verlet is not symplectic. This * Unlike a fixed step size Verlet integrator, variable step size Verlet is not symplectic. This
...@@ -64,20 +64,20 @@ public: ...@@ -64,20 +64,20 @@ public:
* Create a VariableVerletIntegrator. * Create a VariableVerletIntegrator.
* *
* @param stepSize the initial step size to use (in picoseconds) * @param stepSize the initial step size to use (in picoseconds)
* @param accuracy the required accuracy * @param tol the error tolerance
*/ */
VariableVerletIntegrator(double stepSize, double accuracy); VariableVerletIntegrator(double stepSize, double errorTol);
/** /**
* Get the required accuracy. * Get the error tolerance.
*/ */
double getAccuracy() const { double getErrorTolerance() const {
return accuracy; return errorTol;
} }
/** /**
* Set the required accuracy. * Set the error tolerance.
*/ */
void setAccuracy(double acc) { void setErrorTolerance(double tol) {
accuracy = acc; errorTol = tol;
} }
/** /**
* Advance a simulation through time by taking a series of time steps. * Advance a simulation through time by taking a series of time steps.
...@@ -97,7 +97,7 @@ protected: ...@@ -97,7 +97,7 @@ protected:
*/ */
std::vector<std::string> getKernelNames(); std::vector<std::string> getKernelNames();
private: private:
double accuracy; double errorTol;
OpenMMContextImpl* context; OpenMMContextImpl* context;
Kernel kernel; Kernel kernel;
}; };
......
...@@ -39,7 +39,7 @@ using namespace OpenMM; ...@@ -39,7 +39,7 @@ using namespace OpenMM;
using std::string; using std::string;
using std::vector; using std::vector;
VariableVerletIntegrator::VariableVerletIntegrator(double stepSize, double accuracy) : accuracy(accuracy) { VariableVerletIntegrator::VariableVerletIntegrator(double stepSize, double errorTol) : errorTol(errorTol) {
setStepSize(stepSize); setStepSize(stepSize);
setConstraintTolerance(1e-4); setConstraintTolerance(1e-4);
} }
......
...@@ -803,24 +803,24 @@ void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system ...@@ -803,24 +803,24 @@ void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system
void ReferenceIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator) { void ReferenceIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator) {
double stepSize = integrator.getStepSize(); double stepSize = integrator.getStepSize();
double accuracy = integrator.getAccuracy(); double errorTol = integrator.getErrorTolerance();
RealOpenMM** posData = ((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData(); RealOpenMM** posData = ((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData();
RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData(); RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData()); // Reference code needs to be made const correct RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData()); // Reference code needs to be made const correct
if (dynamics == 0 || stepSize != prevStepSize || accuracy != prevAccuracy) { if (dynamics == 0 || stepSize != prevStepSize || errorTol != prevErrorTol) {
// Recreate the computation objects with the new parameters. // Recreate the computation objects with the new parameters.
if (dynamics) { if (dynamics) {
delete dynamics; delete dynamics;
delete constraints; delete constraints;
} }
dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize), accuracy); dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize), errorTol);
vector<ReferenceCCMAAlgorithm::AngleInfo> angles; vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
findAnglesForCCMA(context.getSystem(), angles); findAnglesForCCMA(context.getSystem(), angles);
constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance()); constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
dynamics->setReferenceConstraintAlgorithm(constraints); dynamics->setReferenceConstraintAlgorithm(constraints);
prevStepSize = stepSize; prevStepSize = stepSize;
prevAccuracy = accuracy; prevErrorTol = errorTol;
} }
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses); dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
data.time += dynamics->getLastStepSize(); data.time += dynamics->getLastStepSize();
......
...@@ -472,7 +472,7 @@ private: ...@@ -472,7 +472,7 @@ private:
RealOpenMM* constraintDistances; RealOpenMM* constraintDistances;
int** constraintIndices; int** constraintIndices;
int numConstraints; int numConstraints;
double prevStepSize, prevAccuracy; double prevStepSize, prevErrorTol;
}; };
/** /**
......
...@@ -167,7 +167,6 @@ int ReferenceVariableVerletDynamics::update( int numberOfAtoms, RealOpenMM** ato ...@@ -167,7 +167,6 @@ int ReferenceVariableVerletDynamics::update( int numberOfAtoms, RealOpenMM** ato
// get work arrays // get work arrays
RealOpenMM** xPrime = get2DArrayAtIndex( xPrime2D ); RealOpenMM** xPrime = get2DArrayAtIndex( xPrime2D );
RealOpenMM** vPrime = get2DArrayAtIndex( vPrime2D );
RealOpenMM* inverseMasses = get1DArrayAtIndex( InverseMasses ); RealOpenMM* inverseMasses = get1DArrayAtIndex( InverseMasses );
// first-time-through initialization // first-time-through initialization
...@@ -208,8 +207,8 @@ int ReferenceVariableVerletDynamics::update( int numberOfAtoms, RealOpenMM** ato ...@@ -208,8 +207,8 @@ int ReferenceVariableVerletDynamics::update( int numberOfAtoms, RealOpenMM** ato
for (int i = 0; i < numberOfAtoms; ++i) { for (int i = 0; i < numberOfAtoms; ++i) {
for (int j = 0; j < 3; ++j) { for (int j = 0; j < 3; ++j) {
RealOpenMM xref = atomCoordinates[i][j] + velocities[i][j]*getDeltaT(); RealOpenMM xref = atomCoordinates[i][j] + velocities[i][j]*getDeltaT();
vPrime[i][j] = velocities[i][j] + inverseMasses[i]*forces[i][j]*getDeltaT(); RealOpenMM vPrime = velocities[i][j] + inverseMasses[i]*forces[i][j]*getDeltaT();
xPrime[i][j] = atomCoordinates[i][j] + vPrime[i][j]*getDeltaT(); xPrime[i][j] = atomCoordinates[i][j] + vPrime*getDeltaT();
RealOpenMM xerror = xPrime[i][j]-xref; RealOpenMM xerror = xPrime[i][j]-xref;
error += xerror*xerror; error += xerror*xerror;
} }
......
...@@ -33,7 +33,7 @@ class ReferenceVariableVerletDynamics : public ReferenceDynamics { ...@@ -33,7 +33,7 @@ class ReferenceVariableVerletDynamics : public ReferenceDynamics {
private: private:
enum TwoDArrayIndicies { xPrime2D, vPrime2D, Max2DArrays }; enum TwoDArrayIndicies { xPrime2D, Max2DArrays };
enum OneDArrayIndicies { InverseMasses, Max1DArrays }; enum OneDArrayIndicies { InverseMasses, Max1DArrays };
RealOpenMM _accuracy, _lastStepSize; RealOpenMM _accuracy, _lastStepSize;
......
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