Commit 4d25f4e7 authored by peastman's avatar peastman
Browse files

Major simplification to constraint handling in reference platform

parent 99f3a8a2
...@@ -79,11 +79,12 @@ class ReferenceBrownianDynamics : public ReferenceDynamics { ...@@ -79,11 +79,12 @@ class ReferenceBrownianDynamics : public ReferenceDynamics {
@param velocities velocities @param velocities velocities
@param forces forces @param forces forces
@param masses atom masses @param masses atom masses
@param tolerance the constraint tolerance
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates, void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses ); std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses, RealOpenMM tolerance);
}; };
......
...@@ -37,7 +37,6 @@ class OPENMM_EXPORT ReferenceCCMAAlgorithm : public ReferenceConstraintAlgorithm ...@@ -37,7 +37,6 @@ class OPENMM_EXPORT ReferenceCCMAAlgorithm : public ReferenceConstraintAlgorithm
protected: protected:
int _maximumNumberOfIterations; int _maximumNumberOfIterations;
RealOpenMM _tolerance;
int _numberOfConstraints; int _numberOfConstraints;
std::vector<std::pair<int, int> > _atomIndices; std::vector<std::pair<int, int> > _atomIndices;
...@@ -53,7 +52,7 @@ protected: ...@@ -53,7 +52,7 @@ protected:
private: private:
void applyConstraints(std::vector<OpenMM::RealVec>& atomCoordinates, void applyConstraints(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses, bool constrainingVelocities); std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses, bool constrainingVelocities, RealOpenMM tolerance);
public: public:
class AngleInfo; class AngleInfo;
...@@ -67,9 +66,8 @@ public: ...@@ -67,9 +66,8 @@ public:
* @param distance distances for constraints * @param distance distances for constraints
* @param masses atom masses * @param masses atom masses
* @param angles angle force field terms * @param angles angle force field terms
* @param tolerance constraint tolerance
*/ */
ReferenceCCMAAlgorithm( int numberOfAtoms, int numberOfConstraints, const std::vector<std::pair<int, int> >& atomIndices, const std::vector<RealOpenMM>& distance, std::vector<RealOpenMM>& masses, std::vector<AngleInfo>& angles, RealOpenMM tolerance ); ReferenceCCMAAlgorithm(int numberOfAtoms, int numberOfConstraints, const std::vector<std::pair<int, int> >& atomIndices, const std::vector<RealOpenMM>& distance, std::vector<RealOpenMM>& masses, std::vector<AngleInfo>& angles);
~ReferenceCCMAAlgorithm( ); ~ReferenceCCMAAlgorithm( );
...@@ -88,25 +86,16 @@ public: ...@@ -88,25 +86,16 @@ public:
*/ */
void setMaximumNumberOfIterations( int maximumNumberOfIterations ); void setMaximumNumberOfIterations( int maximumNumberOfIterations );
/**
* Get the constraint tolerance.
*/
RealOpenMM getTolerance( void ) const;
/**
* Set the constraint tolerance.
*/
void setTolerance( RealOpenMM tolerance );
/** /**
* Apply the constraint algorithm. * Apply the constraint algorithm.
* *
* @param atomCoordinates the original atom coordinates * @param atomCoordinates the original atom coordinates
* @param atomCoordinatesP the new atom coordinates * @param atomCoordinatesP the new atom coordinates
* @param inverseMasses 1/mass * @param inverseMasses 1/mass
* @param tolerance the constraint tolerance
*/ */
void apply(std::vector<OpenMM::RealVec>& atomCoordinates, void apply(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses); std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance);
/** /**
* Apply the constraint algorithm to velocities. * Apply the constraint algorithm to velocities.
...@@ -114,9 +103,10 @@ public: ...@@ -114,9 +103,10 @@ public:
* @param atomCoordinates the atom coordinates * @param atomCoordinates the atom coordinates
* @param atomCoordinatesP the velocities to modify * @param atomCoordinatesP the velocities to modify
* @param inverseMasses 1/mass * @param inverseMasses 1/mass
* @param tolerance the constraint tolerance
*/ */
void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates, void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses); std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance);
}; };
class ReferenceCCMAAlgorithm::AngleInfo class ReferenceCCMAAlgorithm::AngleInfo
......
...@@ -36,25 +36,16 @@ public: ...@@ -36,25 +36,16 @@ public:
virtual ~ReferenceConstraintAlgorithm() {}; virtual ~ReferenceConstraintAlgorithm() {};
/**
* Get the constraint tolerance.
*/
virtual RealOpenMM getTolerance() const = 0;
/**
* Set the constraint tolerance.
*/
virtual void setTolerance(RealOpenMM tolerance) = 0;
/** /**
* Apply the constraint algorithm. * Apply the constraint algorithm.
* *
* @param atomCoordinates the original atom coordinates * @param atomCoordinates the original atom coordinates
* @param atomCoordinatesP the new atom coordinates * @param atomCoordinatesP the new atom coordinates
* @param inverseMasses 1/mass * @param inverseMasses 1/mass
* @param tolerance the constraint tolerance
*/ */
virtual void apply(std::vector<OpenMM::RealVec>& atomCoordinates, virtual void apply(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses) = 0; std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) = 0;
/** /**
* Apply the constraint algorithm to velocities. * Apply the constraint algorithm to velocities.
...@@ -62,9 +53,10 @@ public: ...@@ -62,9 +53,10 @@ public:
* @param atomCoordinates the atom coordinates * @param atomCoordinates the atom coordinates
* @param atomCoordinatesP the velocities to modify * @param atomCoordinatesP the velocities to modify
* @param inverseMasses 1/mass * @param inverseMasses 1/mass
* @param tolerance the constraint tolerance
*/ */
virtual void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates, virtual void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses) = 0; std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) = 0;
}; };
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
......
...@@ -46,17 +46,8 @@ namespace OpenMM { ...@@ -46,17 +46,8 @@ namespace OpenMM {
*/ */
class OPENMM_EXPORT ReferenceConstraints : public ReferenceConstraintAlgorithm { class OPENMM_EXPORT ReferenceConstraints : public ReferenceConstraintAlgorithm {
public: public:
ReferenceConstraints(const System& system, RealOpenMM tolerance); ReferenceConstraints(const System& system);
~ReferenceConstraints(); ~ReferenceConstraints();
/**
* Get the constraint tolerance.
*/
RealOpenMM getTolerance() const;
/**
* Set the constraint tolerance.
*/
void setTolerance(RealOpenMM tolerance);
/** /**
* Apply the constraint algorithm. * Apply the constraint algorithm.
...@@ -64,8 +55,9 @@ public: ...@@ -64,8 +55,9 @@ public:
* @param atomCoordinates the original atom coordinates * @param atomCoordinates the original atom coordinates
* @param atomCoordinatesP the new atom coordinates * @param atomCoordinatesP the new atom coordinates
* @param inverseMasses 1/mass * @param inverseMasses 1/mass
* @param tolerance the constraint tolerance
*/ */
void apply(std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses); void apply(std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance);
/** /**
* Apply the constraint algorithm to velocities. * Apply the constraint algorithm to velocities.
...@@ -73,10 +65,10 @@ public: ...@@ -73,10 +65,10 @@ public:
* @param atomCoordinates the atom coordinates * @param atomCoordinates the atom coordinates
* @param atomCoordinatesP the velocities to modify * @param atomCoordinatesP the velocities to modify
* @param inverseMasses 1/mass * @param inverseMasses 1/mass
* @param tolerance the constraint tolerance
*/ */
void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses); void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance);
private: private:
RealOpenMM tolerance;
ReferenceCCMAAlgorithm* ccma; ReferenceCCMAAlgorithm* ccma;
ReferenceSETTLEAlgorithm* settle; ReferenceSETTLEAlgorithm* settle;
}; };
......
...@@ -92,12 +92,13 @@ public: ...@@ -92,12 +92,13 @@ public:
@param globals a map containing values of global variables @param globals a map containing values of global variables
@param perDof the values of per-DOF variables @param perDof the values of per-DOF variables
@param forcesAreValid whether the current forces are valid or need to be recomputed @param forcesAreValid whether the current forces are valid or need to be recomputed
@param tolerance the constraint tolerance
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void update(OpenMM::ContextImpl& context, int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, 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::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); std::map<std::string, RealOpenMM>& globals, std::vector<std::vector<OpenMM::RealVec> >& perDof, bool& forcesAreValid, RealOpenMM tolerance);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
......
...@@ -145,11 +145,12 @@ class OPENMM_EXPORT ReferenceDynamics { ...@@ -145,11 +145,12 @@ class OPENMM_EXPORT ReferenceDynamics {
@param velocities velocities @param velocities velocities
@param forces forces @param forces forces
@param masses atom masses @param masses atom masses
@param tolerance the constraint tolerance
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
virtual void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates, virtual void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses); std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses, RealOpenMM tolerance);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
......
...@@ -191,7 +191,7 @@ private: ...@@ -191,7 +191,7 @@ private:
class ReferenceApplyConstraintsKernel : public ApplyConstraintsKernel { class ReferenceApplyConstraintsKernel : public ApplyConstraintsKernel {
public: public:
ReferenceApplyConstraintsKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : ReferenceApplyConstraintsKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) :
ApplyConstraintsKernel(name, platform), data(data), constraints(0) { ApplyConstraintsKernel(name, platform), data(data) {
} }
~ReferenceApplyConstraintsKernel(); ~ReferenceApplyConstraintsKernel();
/** /**
...@@ -216,12 +216,8 @@ public: ...@@ -216,12 +216,8 @@ public:
void applyToVelocities(ContextImpl& context, double tol); void applyToVelocities(ContextImpl& context, double tol);
private: private:
ReferencePlatform::PlatformData& data; ReferencePlatform::PlatformData& data;
ReferenceConstraintAlgorithm* constraints;
std::vector<RealOpenMM> masses; std::vector<RealOpenMM> masses;
std::vector<RealOpenMM> inverseMasses; std::vector<RealOpenMM> inverseMasses;
std::vector<std::pair<int, int> > constraintIndices;
std::vector<RealOpenMM> constraintDistances;
int numConstraints;
}; };
/** /**
...@@ -871,7 +867,7 @@ private: ...@@ -871,7 +867,7 @@ private:
class ReferenceIntegrateVerletStepKernel : public IntegrateVerletStepKernel { class ReferenceIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
public: public:
ReferenceIntegrateVerletStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVerletStepKernel(name, platform), ReferenceIntegrateVerletStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVerletStepKernel(name, platform),
data(data), dynamics(0), constraints(0) { data(data), dynamics(0) {
} }
~ReferenceIntegrateVerletStepKernel(); ~ReferenceIntegrateVerletStepKernel();
/** /**
...@@ -898,9 +894,7 @@ public: ...@@ -898,9 +894,7 @@ public:
private: private:
ReferencePlatform::PlatformData& data; ReferencePlatform::PlatformData& data;
ReferenceVerletDynamics* dynamics; ReferenceVerletDynamics* dynamics;
ReferenceConstraintAlgorithm* constraints;
std::vector<RealOpenMM> masses; std::vector<RealOpenMM> masses;
int numConstraints;
double prevStepSize; double prevStepSize;
}; };
...@@ -910,7 +904,7 @@ private: ...@@ -910,7 +904,7 @@ private:
class ReferenceIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel { class ReferenceIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
public: public:
ReferenceIntegrateLangevinStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateLangevinStepKernel(name, platform), ReferenceIntegrateLangevinStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateLangevinStepKernel(name, platform),
data(data), dynamics(0), constraints(0) { data(data), dynamics(0) {
} }
~ReferenceIntegrateLangevinStepKernel(); ~ReferenceIntegrateLangevinStepKernel();
/** /**
...@@ -937,9 +931,7 @@ public: ...@@ -937,9 +931,7 @@ public:
private: private:
ReferencePlatform::PlatformData& data; ReferencePlatform::PlatformData& data;
ReferenceStochasticDynamics* dynamics; ReferenceStochasticDynamics* dynamics;
ReferenceConstraintAlgorithm* constraints;
std::vector<RealOpenMM> masses; std::vector<RealOpenMM> masses;
int numConstraints;
double prevTemp, prevFriction, prevStepSize; double prevTemp, prevFriction, prevStepSize;
}; };
...@@ -949,7 +941,7 @@ private: ...@@ -949,7 +941,7 @@ private:
class ReferenceIntegrateBrownianStepKernel : public IntegrateBrownianStepKernel { class ReferenceIntegrateBrownianStepKernel : public IntegrateBrownianStepKernel {
public: public:
ReferenceIntegrateBrownianStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateBrownianStepKernel(name, platform), ReferenceIntegrateBrownianStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateBrownianStepKernel(name, platform),
data(data), dynamics(0), constraints(0) { data(data), dynamics(0) {
} }
~ReferenceIntegrateBrownianStepKernel(); ~ReferenceIntegrateBrownianStepKernel();
/** /**
...@@ -976,9 +968,7 @@ public: ...@@ -976,9 +968,7 @@ public:
private: private:
ReferencePlatform::PlatformData& data; ReferencePlatform::PlatformData& data;
ReferenceBrownianDynamics* dynamics; ReferenceBrownianDynamics* dynamics;
ReferenceConstraintAlgorithm* constraints;
std::vector<RealOpenMM> masses; std::vector<RealOpenMM> masses;
int numConstraints;
double prevTemp, prevFriction, prevStepSize; double prevTemp, prevFriction, prevStepSize;
}; };
...@@ -988,7 +978,7 @@ private: ...@@ -988,7 +978,7 @@ private:
class ReferenceIntegrateVariableLangevinStepKernel : public IntegrateVariableLangevinStepKernel { class ReferenceIntegrateVariableLangevinStepKernel : public IntegrateVariableLangevinStepKernel {
public: public:
ReferenceIntegrateVariableLangevinStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVariableLangevinStepKernel(name, platform), ReferenceIntegrateVariableLangevinStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVariableLangevinStepKernel(name, platform),
data(data), dynamics(0), constraints(0) { data(data), dynamics(0) {
} }
~ReferenceIntegrateVariableLangevinStepKernel(); ~ReferenceIntegrateVariableLangevinStepKernel();
/** /**
...@@ -1017,9 +1007,7 @@ public: ...@@ -1017,9 +1007,7 @@ public:
private: private:
ReferencePlatform::PlatformData& data; ReferencePlatform::PlatformData& data;
ReferenceVariableStochasticDynamics* dynamics; ReferenceVariableStochasticDynamics* dynamics;
ReferenceConstraintAlgorithm* constraints;
std::vector<RealOpenMM> masses; std::vector<RealOpenMM> masses;
int numConstraints;
double prevTemp, prevFriction, prevErrorTol; double prevTemp, prevFriction, prevErrorTol;
}; };
...@@ -1029,7 +1017,7 @@ private: ...@@ -1029,7 +1017,7 @@ private:
class ReferenceIntegrateVariableVerletStepKernel : public IntegrateVariableVerletStepKernel { class ReferenceIntegrateVariableVerletStepKernel : public IntegrateVariableVerletStepKernel {
public: public:
ReferenceIntegrateVariableVerletStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVariableVerletStepKernel(name, platform), ReferenceIntegrateVariableVerletStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVariableVerletStepKernel(name, platform),
data(data), dynamics(0), constraints(0) { data(data), dynamics(0) {
} }
~ReferenceIntegrateVariableVerletStepKernel(); ~ReferenceIntegrateVariableVerletStepKernel();
/** /**
...@@ -1058,9 +1046,7 @@ public: ...@@ -1058,9 +1046,7 @@ public:
private: private:
ReferencePlatform::PlatformData& data; ReferencePlatform::PlatformData& data;
ReferenceVariableVerletDynamics* dynamics; ReferenceVariableVerletDynamics* dynamics;
ReferenceConstraintAlgorithm* constraints;
std::vector<RealOpenMM> masses; std::vector<RealOpenMM> masses;
int numConstraints;
double prevErrorTol; double prevErrorTol;
}; };
...@@ -1070,7 +1056,7 @@ private: ...@@ -1070,7 +1056,7 @@ private:
class ReferenceIntegrateCustomStepKernel : public IntegrateCustomStepKernel { class ReferenceIntegrateCustomStepKernel : public IntegrateCustomStepKernel {
public: public:
ReferenceIntegrateCustomStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateCustomStepKernel(name, platform), ReferenceIntegrateCustomStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateCustomStepKernel(name, platform),
data(data), dynamics(0), constraints(0) { data(data), dynamics(0) {
} }
~ReferenceIntegrateCustomStepKernel(); ~ReferenceIntegrateCustomStepKernel();
/** /**
...@@ -1135,10 +1121,8 @@ public: ...@@ -1135,10 +1121,8 @@ public:
private: private:
ReferencePlatform::PlatformData& data; ReferencePlatform::PlatformData& data;
ReferenceCustomDynamics* dynamics; ReferenceCustomDynamics* dynamics;
ReferenceConstraintAlgorithm* constraints;
std::vector<RealOpenMM> masses, globalValues; std::vector<RealOpenMM> masses, globalValues;
std::vector<std::vector<OpenMM::RealVec> > perDofValues; std::vector<std::vector<OpenMM::RealVec> > perDofValues;
int numConstraints;
}; };
/** /**
......
...@@ -33,6 +33,7 @@ ...@@ -33,6 +33,7 @@
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
#include "openmm/Platform.h" #include "openmm/Platform.h"
#include "openmm/System.h"
#include "openmm/internal/windowsExport.h" #include "openmm/internal/windowsExport.h"
namespace OpenMM { namespace OpenMM {
...@@ -57,7 +58,7 @@ public: ...@@ -57,7 +58,7 @@ public:
class ReferencePlatform::PlatformData { class ReferencePlatform::PlatformData {
public: public:
PlatformData(int numParticles); PlatformData(const System& system);
~PlatformData(); ~PlatformData();
int numParticles, stepCount; int numParticles, stepCount;
double time; double time;
...@@ -65,6 +66,7 @@ public: ...@@ -65,6 +66,7 @@ public:
void* velocities; void* velocities;
void* forces; void* forces;
void* periodicBoxSize; void* periodicBoxSize;
void* constraints;
}; };
} // namespace OpenMM } // namespace OpenMM
......
...@@ -44,16 +44,7 @@ namespace OpenMM { ...@@ -44,16 +44,7 @@ namespace OpenMM {
class OPENMM_EXPORT ReferenceSETTLEAlgorithm : public ReferenceConstraintAlgorithm { class OPENMM_EXPORT ReferenceSETTLEAlgorithm : public ReferenceConstraintAlgorithm {
public: public:
ReferenceSETTLEAlgorithm(const std::vector<int>& atom1, const std::vector<int>& atom2, const std::vector<int>& atom3, ReferenceSETTLEAlgorithm(const std::vector<int>& atom1, const std::vector<int>& atom2, const std::vector<int>& atom3,
const std::vector<RealOpenMM>& distance1, const std::vector<RealOpenMM>& distance2, std::vector<RealOpenMM>& masses, RealOpenMM tolerance); const std::vector<RealOpenMM>& distance1, const std::vector<RealOpenMM>& distance2, std::vector<RealOpenMM>& masses);
/**
* Get the constraint tolerance.
*/
RealOpenMM getTolerance() const;
/**
* Set the constraint tolerance.
*/
void setTolerance(RealOpenMM tolerance);
/** /**
* Apply the constraint algorithm. * Apply the constraint algorithm.
...@@ -61,8 +52,9 @@ public: ...@@ -61,8 +52,9 @@ public:
* @param atomCoordinates the original atom coordinates * @param atomCoordinates the original atom coordinates
* @param atomCoordinatesP the new atom coordinates * @param atomCoordinatesP the new atom coordinates
* @param inverseMasses 1/mass * @param inverseMasses 1/mass
* @param tolerance the constraint tolerance
*/ */
void apply(std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses); void apply(std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance);
/** /**
* Apply the constraint algorithm to velocities. * Apply the constraint algorithm to velocities.
...@@ -70,10 +62,10 @@ public: ...@@ -70,10 +62,10 @@ public:
* @param atomCoordinates the atom coordinates * @param atomCoordinates the atom coordinates
* @param atomCoordinatesP the velocities to modify * @param atomCoordinatesP the velocities to modify
* @param inverseMasses 1/mass * @param inverseMasses 1/mass
* @param tolerance the constraint tolerance
*/ */
void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses); void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance);
private: private:
RealOpenMM tolerance;
std::vector<int> atom1; std::vector<int> atom1;
std::vector<int> atom2; std::vector<int> atom2;
std::vector<int> atom3; std::vector<int> atom3;
......
...@@ -79,11 +79,12 @@ class ReferenceStochasticDynamics : public ReferenceDynamics { ...@@ -79,11 +79,12 @@ class ReferenceStochasticDynamics : public ReferenceDynamics {
@param velocities velocities @param velocities velocities
@param forces forces @param forces forces
@param masses atom masses @param masses atom masses
@param tolerance the constraint tolerance
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates, void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses ); std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses, RealOpenMM tolerance);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
......
...@@ -98,11 +98,12 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics { ...@@ -98,11 +98,12 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
@param forces forces @param forces forces
@param masses atom masses @param masses atom masses
@param maxStepSize maximum time step @param maxStepSize maximum time step
@param tolerance the constraint tolerance
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates, void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses, RealOpenMM maxStepSize ); std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses, RealOpenMM maxStepSize, RealOpenMM tolerance);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
......
...@@ -86,11 +86,12 @@ class ReferenceVariableVerletDynamics : public ReferenceDynamics { ...@@ -86,11 +86,12 @@ class ReferenceVariableVerletDynamics : public ReferenceDynamics {
@param forces forces @param forces forces
@param masses atom masses @param masses atom masses
@param maxStepSize maximum time step @param maxStepSize maximum time step
@param tolerance the constraint tolerance
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates, void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses, RealOpenMM maxStepSize ); std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses, RealOpenMM maxStepSize, RealOpenMM tolerance);
}; };
......
...@@ -68,11 +68,12 @@ class ReferenceVerletDynamics : public ReferenceDynamics { ...@@ -68,11 +68,12 @@ class ReferenceVerletDynamics : public ReferenceDynamics {
@param velocities velocities @param velocities velocities
@param forces forces @param forces forces
@param masses atom masses @param masses atom masses
@param tolerance the constraint tolerance
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates, void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses ); std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses, RealOpenMM tolerance);
}; };
......
...@@ -133,11 +133,16 @@ static RealVec& extractBoxSize(ContextImpl& context) { ...@@ -133,11 +133,16 @@ static RealVec& extractBoxSize(ContextImpl& context) {
return *(RealVec*) data->periodicBoxSize; return *(RealVec*) data->periodicBoxSize;
} }
static ReferenceConstraints& extractConstraints(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *(ReferenceConstraints*) data->constraints;
}
/** /**
* Compute the kinetic energy of the system, possibly shifting the velocities in time to account * Compute the kinetic energy of the system, possibly shifting the velocities in time to account
* for a leapfrog integrator. * for a leapfrog integrator.
*/ */
static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& masses, double timeShift, ReferenceConstraintAlgorithm* constraints) { static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& masses, double timeShift) {
vector<RealVec>& posData = extractPositions(context); vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& velData = extractVelocities(context); vector<RealVec>& velData = extractVelocities(context);
vector<RealVec>& forceData = extractForces(context); vector<RealVec>& forceData = extractForces(context);
...@@ -155,13 +160,10 @@ static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& ...@@ -155,13 +160,10 @@ static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>&
// Apply constraints to them. // Apply constraints to them.
if (constraints != NULL) {
vector<double> inverseMasses(numParticles); vector<double> inverseMasses(numParticles);
for (int i = 0; i < numParticles; i++) for (int i = 0; i < numParticles; i++)
inverseMasses[i] = (masses[i] == 0 ? 0 : 1/masses[i]); inverseMasses[i] = (masses[i] == 0 ? 0 : 1/masses[i]);
constraints->setTolerance(1e-4); extractConstraints(context).applyToVelocities(posData, shiftedVel, inverseMasses, 1e-4);
constraints->applyToVelocities(posData, shiftedVel, inverseMasses);
}
// Compute the kinetic energy. // Compute the kinetic energy.
...@@ -302,39 +304,21 @@ void ReferenceApplyConstraintsKernel::initialize(const System& system) { ...@@ -302,39 +304,21 @@ void ReferenceApplyConstraintsKernel::initialize(const System& system) {
masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i)); masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
inverseMasses[i] = 1.0/masses[i]; inverseMasses[i] = 1.0/masses[i];
} }
for (int i = 0; i < system.getNumConstraints(); ++i) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
constraintIndices.push_back(make_pair(particle1, particle2));
constraintDistances.push_back(distance);
}
}
numConstraints = constraintIndices.size();
} }
ReferenceApplyConstraintsKernel::~ReferenceApplyConstraintsKernel() { ReferenceApplyConstraintsKernel::~ReferenceApplyConstraintsKernel() {
if (constraints)
delete constraints;
} }
void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) { void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
if (constraints == NULL)
constraints = new ReferenceConstraints(context.getSystem(), (RealOpenMM) tol);
vector<RealVec>& positions = extractPositions(context); vector<RealVec>& positions = extractPositions(context);
constraints->setTolerance(tol); extractConstraints(context).apply(positions, positions, inverseMasses, tol);
constraints->apply(positions, positions, inverseMasses);
ReferenceVirtualSites::computePositions(context.getSystem(), positions); ReferenceVirtualSites::computePositions(context.getSystem(), positions);
} }
void ReferenceApplyConstraintsKernel::applyToVelocities(ContextImpl& context, double tol) { void ReferenceApplyConstraintsKernel::applyToVelocities(ContextImpl& context, double tol) {
if (constraints == NULL)
constraints = new ReferenceConstraints(context.getSystem(), (RealOpenMM) tol);
vector<RealVec>& positions = extractPositions(context); vector<RealVec>& positions = extractPositions(context);
vector<RealVec>& velocities = extractVelocities(context); vector<RealVec>& velocities = extractVelocities(context);
constraints->setTolerance(tol); extractConstraints(context).applyToVelocities(positions, velocities, inverseMasses, tol);
constraints->applyToVelocities(positions, velocities, inverseMasses);
} }
void ReferenceVirtualSitesKernel::initialize(const System& system) { void ReferenceVirtualSitesKernel::initialize(const System& system) {
...@@ -1686,8 +1670,6 @@ void ReferenceCalcCustomCompoundBondForceKernel::copyParametersToContext(Context ...@@ -1686,8 +1670,6 @@ void ReferenceCalcCustomCompoundBondForceKernel::copyParametersToContext(Context
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() { ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
if (dynamics) if (dynamics)
delete dynamics; delete dynamics;
if (constraints)
delete constraints;
} }
void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) { void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
...@@ -1695,19 +1677,6 @@ void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const ...@@ -1695,19 +1677,6 @@ void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const
masses.resize(numParticles); masses.resize(numParticles);
for (int i = 0; i < numParticles; ++i) for (int i = 0; i < numParticles; ++i)
masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i)); masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
vector<pair<int, int> > constraintIndices;
vector<RealOpenMM> constraintDistances;
for (int i = 0; i < system.getNumConstraints(); ++i) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
constraintIndices.push_back(make_pair(particle1, particle2));
constraintDistances.push_back(distance);
}
}
numConstraints = constraintIndices.size();
constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
} }
void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) { void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
...@@ -1721,24 +1690,21 @@ void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const Ver ...@@ -1721,24 +1690,21 @@ void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const Ver
if (dynamics) if (dynamics)
delete dynamics; delete dynamics;
dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize) ); dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize) );
dynamics->setReferenceConstraintAlgorithm(constraints); dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevStepSize = stepSize; prevStepSize = stepSize;
} }
constraints->setTolerance(integrator.getConstraintTolerance()); dynamics->update(context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
dynamics->update(context.getSystem(), posData, velData, forceData, masses);
data.time += stepSize; data.time += stepSize;
data.stepCount++; data.stepCount++;
} }
double ReferenceIntegrateVerletStepKernel::computeKineticEnergy(ContextImpl& context, const VerletIntegrator& integrator) { double ReferenceIntegrateVerletStepKernel::computeKineticEnergy(ContextImpl& context, const VerletIntegrator& integrator) {
return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize(), constraints); return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
} }
ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() { ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
if (dynamics) if (dynamics)
delete dynamics; delete dynamics;
if (constraints)
delete constraints;
} }
void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) { void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
...@@ -1746,20 +1712,7 @@ void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, cons ...@@ -1746,20 +1712,7 @@ void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, cons
masses.resize(numParticles); masses.resize(numParticles);
for (int i = 0; i < numParticles; ++i) for (int i = 0; i < numParticles; ++i)
masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i)); masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
vector<pair<int, int> > constraintIndices;
vector<RealOpenMM> constraintDistances;
for (int i = 0; i < system.getNumConstraints(); ++i) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
constraintIndices.push_back(make_pair(particle1, particle2));
constraintDistances.push_back(distance);
}
}
numConstraints = constraintIndices.size();
SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed()); SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
} }
void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) { void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
...@@ -1780,26 +1733,23 @@ void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const L ...@@ -1780,26 +1733,23 @@ void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const L
static_cast<RealOpenMM>(stepSize), static_cast<RealOpenMM>(stepSize),
static_cast<RealOpenMM>(tau), static_cast<RealOpenMM>(tau),
static_cast<RealOpenMM>(temperature) ); static_cast<RealOpenMM>(temperature) );
dynamics->setReferenceConstraintAlgorithm(constraints); dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevTemp = temperature; prevTemp = temperature;
prevFriction = friction; prevFriction = friction;
prevStepSize = stepSize; prevStepSize = stepSize;
} }
constraints->setTolerance(integrator.getConstraintTolerance()); dynamics->update(context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
dynamics->update(context.getSystem(), posData, velData, forceData, masses);
data.time += stepSize; data.time += stepSize;
data.stepCount++; data.stepCount++;
} }
double ReferenceIntegrateLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const LangevinIntegrator& integrator) { double ReferenceIntegrateLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const LangevinIntegrator& integrator) {
return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize(), constraints); return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
} }
ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() { ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
if (dynamics) if (dynamics)
delete dynamics; delete dynamics;
if (constraints)
delete constraints;
} }
void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) { void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
...@@ -1807,20 +1757,7 @@ void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, cons ...@@ -1807,20 +1757,7 @@ void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, cons
masses.resize(numParticles); masses.resize(numParticles);
for (int i = 0; i < numParticles; ++i) for (int i = 0; i < numParticles; ++i)
masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i)); masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
vector<pair<int, int> > constraintIndices;
vector<RealOpenMM> constraintDistances;
for (int i = 0; i < system.getNumConstraints(); ++i) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
constraintIndices.push_back(make_pair(particle1, particle2));
constraintDistances.push_back(distance);
}
}
numConstraints = constraintIndices.size();
SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed()); SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
} }
void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) { void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
...@@ -1840,26 +1777,23 @@ void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const B ...@@ -1840,26 +1777,23 @@ void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const B
static_cast<RealOpenMM>(stepSize), static_cast<RealOpenMM>(stepSize),
static_cast<RealOpenMM>(friction), static_cast<RealOpenMM>(friction),
static_cast<RealOpenMM>(temperature) ); static_cast<RealOpenMM>(temperature) );
dynamics->setReferenceConstraintAlgorithm(constraints); dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevTemp = temperature; prevTemp = temperature;
prevFriction = friction; prevFriction = friction;
prevStepSize = stepSize; prevStepSize = stepSize;
} }
constraints->setTolerance(integrator.getConstraintTolerance()); dynamics->update(context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
dynamics->update(context.getSystem(), posData, velData, forceData, masses);
data.time += stepSize; data.time += stepSize;
data.stepCount++; data.stepCount++;
} }
double ReferenceIntegrateBrownianStepKernel::computeKineticEnergy(ContextImpl& context, const BrownianIntegrator& integrator) { double ReferenceIntegrateBrownianStepKernel::computeKineticEnergy(ContextImpl& context, const BrownianIntegrator& integrator) {
return computeShiftedKineticEnergy(context, masses, 0, constraints); return computeShiftedKineticEnergy(context, masses, 0);
} }
ReferenceIntegrateVariableLangevinStepKernel::~ReferenceIntegrateVariableLangevinStepKernel() { ReferenceIntegrateVariableLangevinStepKernel::~ReferenceIntegrateVariableLangevinStepKernel() {
if (dynamics) if (dynamics)
delete dynamics; delete dynamics;
if (constraints)
delete constraints;
} }
void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& system, const VariableLangevinIntegrator& integrator) { void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& system, const VariableLangevinIntegrator& integrator) {
...@@ -1867,20 +1801,7 @@ void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& syst ...@@ -1867,20 +1801,7 @@ void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& syst
masses.resize(numParticles); masses.resize(numParticles);
for (int i = 0; i < numParticles; ++i) for (int i = 0; i < numParticles; ++i)
masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i)); masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
vector<pair<int, int> > constraintIndices;
vector<RealOpenMM> constraintDistances;
for (int i = 0; i < system.getNumConstraints(); ++i) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
constraintIndices.push_back(make_pair(particle1, particle2));
constraintDistances.push_back(distance);
}
}
numConstraints = constraintIndices.size();
SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed()); SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
} }
double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) { double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
...@@ -1897,14 +1818,13 @@ double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& contex ...@@ -1897,14 +1818,13 @@ double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& contex
delete dynamics; delete dynamics;
RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction ); RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
dynamics = new ReferenceVariableStochasticDynamics(context.getSystem().getNumParticles(), (RealOpenMM) tau, (RealOpenMM) temperature, (RealOpenMM) errorTol); dynamics = new ReferenceVariableStochasticDynamics(context.getSystem().getNumParticles(), (RealOpenMM) tau, (RealOpenMM) temperature, (RealOpenMM) errorTol);
dynamics->setReferenceConstraintAlgorithm(constraints); dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevTemp = temperature; prevTemp = temperature;
prevFriction = friction; prevFriction = friction;
prevErrorTol = errorTol; prevErrorTol = errorTol;
} }
constraints->setTolerance(integrator.getConstraintTolerance());
RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time); RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize); dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize, integrator.getConstraintTolerance());
data.time += dynamics->getDeltaT(); data.time += dynamics->getDeltaT();
if (dynamics->getDeltaT() == maxStepSize) if (dynamics->getDeltaT() == maxStepSize)
data.time = maxTime; // Avoid round-off error data.time = maxTime; // Avoid round-off error
...@@ -1913,14 +1833,12 @@ double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& contex ...@@ -1913,14 +1833,12 @@ double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& contex
} }
double ReferenceIntegrateVariableLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const VariableLangevinIntegrator& integrator) { double ReferenceIntegrateVariableLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const VariableLangevinIntegrator& integrator) {
return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize(), constraints); return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
} }
ReferenceIntegrateVariableVerletStepKernel::~ReferenceIntegrateVariableVerletStepKernel() { ReferenceIntegrateVariableVerletStepKernel::~ReferenceIntegrateVariableVerletStepKernel() {
if (dynamics) if (dynamics)
delete dynamics; delete dynamics;
if (constraints)
delete constraints;
} }
void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system, const VariableVerletIntegrator& integrator) { void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system, const VariableVerletIntegrator& integrator) {
...@@ -1928,19 +1846,6 @@ void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system ...@@ -1928,19 +1846,6 @@ void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system
masses.resize(numParticles); masses.resize(numParticles);
for (int i = 0; i < numParticles; ++i) for (int i = 0; i < numParticles; ++i)
masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i)); masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
vector<pair<int, int> > constraintIndices;
vector<RealOpenMM> constraintDistances;
for (int i = 0; i < system.getNumConstraints(); ++i) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
constraintIndices.push_back(make_pair(particle1, particle2));
constraintDistances.push_back(distance);
}
}
numConstraints = constraintIndices.size();
constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
} }
double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) { double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
...@@ -1954,12 +1859,11 @@ double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, ...@@ -1954,12 +1859,11 @@ double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context,
if (dynamics) if (dynamics)
delete dynamics; delete dynamics;
dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), (RealOpenMM) errorTol); dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), (RealOpenMM) errorTol);
dynamics->setReferenceConstraintAlgorithm(constraints); dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevErrorTol = errorTol; prevErrorTol = errorTol;
} }
constraints->setTolerance(integrator.getConstraintTolerance());
RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time); RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize); dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize, integrator.getConstraintTolerance());
data.time += dynamics->getDeltaT(); data.time += dynamics->getDeltaT();
if (dynamics->getDeltaT() == maxStepSize) if (dynamics->getDeltaT() == maxStepSize)
data.time = maxTime; // Avoid round-off error data.time = maxTime; // Avoid round-off error
...@@ -1968,14 +1872,12 @@ double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, ...@@ -1968,14 +1872,12 @@ double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context,
} }
double ReferenceIntegrateVariableVerletStepKernel::computeKineticEnergy(ContextImpl& context, const VariableVerletIntegrator& integrator) { double ReferenceIntegrateVariableVerletStepKernel::computeKineticEnergy(ContextImpl& context, const VariableVerletIntegrator& integrator) {
return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize(), constraints); return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
} }
ReferenceIntegrateCustomStepKernel::~ReferenceIntegrateCustomStepKernel() { ReferenceIntegrateCustomStepKernel::~ReferenceIntegrateCustomStepKernel() {
if (dynamics) if (dynamics)
delete dynamics; delete dynamics;
if (constraints)
delete constraints;
} }
void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const CustomIntegrator& integrator) { void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const CustomIntegrator& integrator) {
...@@ -1983,18 +1885,6 @@ void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const ...@@ -1983,18 +1885,6 @@ void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const
masses.resize(numParticles); masses.resize(numParticles);
for (int i = 0; i < numParticles; ++i) for (int i = 0; i < numParticles; ++i)
masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i)); masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
vector<pair<int, int> > constraintIndices;
vector<RealOpenMM> constraintDistances;
for (int i = 0; i < system.getNumConstraints(); ++i) {
int particle1, particle2;
double distance;
system.getConstraintParameters(i, particle1, particle2, distance);
if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
constraintIndices.push_back(make_pair(particle1, particle2));
constraintDistances.push_back(distance);
}
}
numConstraints = constraintIndices.size();
perDofValues.resize(integrator.getNumPerDofVariables()); perDofValues.resize(integrator.getNumPerDofVariables());
for (int i = 0; i < (int) perDofValues.size(); i++) for (int i = 0; i < (int) perDofValues.size(); i++)
perDofValues[i].resize(numParticles); perDofValues[i].resize(numParticles);
...@@ -2003,8 +1893,6 @@ void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const ...@@ -2003,8 +1893,6 @@ void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const
dynamics = new ReferenceCustomDynamics(system.getNumParticles(), integrator); dynamics = new ReferenceCustomDynamics(system.getNumParticles(), integrator);
SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed()); SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
dynamics->setReferenceConstraintAlgorithm(constraints);
} }
void ReferenceIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid) { void ReferenceIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid) {
...@@ -2021,8 +1909,8 @@ void ReferenceIntegrateCustomStepKernel::execute(ContextImpl& context, CustomInt ...@@ -2021,8 +1909,8 @@ void ReferenceIntegrateCustomStepKernel::execute(ContextImpl& context, CustomInt
// Execute the step. // Execute the step.
constraints->setTolerance(integrator.getConstraintTolerance()); dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
dynamics->update(context, context.getSystem().getNumParticles(), posData, velData, forceData, masses, globals, perDofValues, forcesAreValid); dynamics->update(context, context.getSystem().getNumParticles(), posData, velData, forceData, masses, globals, perDofValues, forcesAreValid, integrator.getConstraintTolerance());
// Record changed global variables. // Record changed global variables.
......
...@@ -30,6 +30,7 @@ ...@@ -30,6 +30,7 @@
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
#include "ReferencePlatform.h" #include "ReferencePlatform.h"
#include "ReferenceConstraints.h"
#include "ReferenceKernelFactory.h" #include "ReferenceKernelFactory.h"
#include "ReferenceKernels.h" #include "ReferenceKernels.h"
#include "openmm/internal/ContextImpl.h" #include "openmm/internal/ContextImpl.h"
...@@ -82,7 +83,7 @@ bool ReferencePlatform::supportsDoublePrecision() const { ...@@ -82,7 +83,7 @@ bool ReferencePlatform::supportsDoublePrecision() const {
} }
void ReferencePlatform::contextCreated(ContextImpl& context, const map<string, string>& properties) const { void ReferencePlatform::contextCreated(ContextImpl& context, const map<string, string>& properties) const {
context.setPlatformData(new PlatformData(context.getSystem().getNumParticles())); context.setPlatformData(new PlatformData(context.getSystem()));
} }
void ReferencePlatform::contextDestroyed(ContextImpl& context) const { void ReferencePlatform::contextDestroyed(ContextImpl& context) const {
...@@ -90,11 +91,12 @@ void ReferencePlatform::contextDestroyed(ContextImpl& context) const { ...@@ -90,11 +91,12 @@ void ReferencePlatform::contextDestroyed(ContextImpl& context) const {
delete data; delete data;
} }
ReferencePlatform::PlatformData::PlatformData(int numParticles) : time(0.0), stepCount(0), numParticles(numParticles) { ReferencePlatform::PlatformData::PlatformData(const System& system) : time(0.0), stepCount(0), numParticles(system.getNumParticles()) {
positions = new vector<RealVec>(numParticles); positions = new vector<RealVec>(numParticles);
velocities = new vector<RealVec>(numParticles); velocities = new vector<RealVec>(numParticles);
forces = new vector<RealVec>(numParticles); forces = new vector<RealVec>(numParticles);
periodicBoxSize = new RealVec(); periodicBoxSize = new RealVec();
constraints = new ReferenceConstraints(system);
} }
ReferencePlatform::PlatformData::~PlatformData() { ReferencePlatform::PlatformData::~PlatformData() {
...@@ -102,4 +104,5 @@ ReferencePlatform::PlatformData::~PlatformData() { ...@@ -102,4 +104,5 @@ ReferencePlatform::PlatformData::~PlatformData() {
delete (vector<RealVec>*) velocities; delete (vector<RealVec>*) velocities;
delete (vector<RealVec>*) forces; delete (vector<RealVec>*) forces;
delete (RealVec*) periodicBoxSize; delete (RealVec*) periodicBoxSize;
delete (ReferenceConstraints*) constraints;
} }
...@@ -125,7 +125,7 @@ RealOpenMM ReferenceBrownianDynamics::getFriction( void ) const { ...@@ -125,7 +125,7 @@ RealOpenMM ReferenceBrownianDynamics::getFriction( void ) const {
void ReferenceBrownianDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates, void ReferenceBrownianDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities, vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& masses ){ vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM tolerance) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -162,7 +162,7 @@ void ReferenceBrownianDynamics::update(const OpenMM::System& system, vector<Real ...@@ -162,7 +162,7 @@ void ReferenceBrownianDynamics::update(const OpenMM::System& system, vector<Real
} }
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm(); ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if (referenceConstraintAlgorithm) if (referenceConstraintAlgorithm)
referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses); referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses, tolerance);
// Update the positions and velocities. // Update the positions and velocities.
......
...@@ -48,14 +48,12 @@ ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm(int numberOfAtoms, ...@@ -48,14 +48,12 @@ ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm(int numberOfAtoms,
const vector<pair<int, int> >& atomIndices, const vector<pair<int, int> >& atomIndices,
const vector<RealOpenMM>& distance, const vector<RealOpenMM>& distance,
vector<RealOpenMM>& masses, vector<RealOpenMM>& masses,
vector<AngleInfo>& angles, vector<AngleInfo>& angles) {
RealOpenMM tolerance) {
_numberOfConstraints = numberOfConstraints; _numberOfConstraints = numberOfConstraints;
_atomIndices = atomIndices; _atomIndices = atomIndices;
_distance = distance; _distance = distance;
_maximumNumberOfIterations = 150; _maximumNumberOfIterations = 150;
_tolerance = tolerance;
_hasInitializedMasses = false; _hasInitializedMasses = false;
// work arrays // work arrays
...@@ -201,28 +199,20 @@ void ReferenceCCMAAlgorithm::setMaximumNumberOfIterations(int maximumNumberOfIte ...@@ -201,28 +199,20 @@ void ReferenceCCMAAlgorithm::setMaximumNumberOfIterations(int maximumNumberOfIte
_maximumNumberOfIterations = maximumNumberOfIterations; _maximumNumberOfIterations = maximumNumberOfIterations;
} }
RealOpenMM ReferenceCCMAAlgorithm::getTolerance() const {
return _tolerance;
}
void ReferenceCCMAAlgorithm::setTolerance(RealOpenMM tolerance) {
_tolerance = tolerance;
}
void ReferenceCCMAAlgorithm::apply(vector<RealVec>& atomCoordinates, void ReferenceCCMAAlgorithm::apply(vector<RealVec>& atomCoordinates,
vector<RealVec>& atomCoordinatesP, vector<RealVec>& atomCoordinatesP,
vector<RealOpenMM>& inverseMasses) { vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) {
applyConstraints(atomCoordinates, atomCoordinatesP, inverseMasses, false); applyConstraints(atomCoordinates, atomCoordinatesP, inverseMasses, false, tolerance);
} }
void ReferenceCCMAAlgorithm::applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates, void ReferenceCCMAAlgorithm::applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses) { std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) {
applyConstraints(atomCoordinates, velocities, inverseMasses, true); applyConstraints(atomCoordinates, velocities, inverseMasses, true, tolerance);
} }
void ReferenceCCMAAlgorithm::applyConstraints(vector<RealVec>& atomCoordinates, void ReferenceCCMAAlgorithm::applyConstraints(vector<RealVec>& atomCoordinates,
vector<RealVec>& atomCoordinatesP, vector<RealVec>& atomCoordinatesP,
vector<RealOpenMM>& inverseMasses, bool constrainingVelocities) { vector<RealOpenMM>& inverseMasses, bool constrainingVelocities, RealOpenMM tolerance) {
// temp arrays // temp arrays
vector<RealVec>& r_ij = _r_ij; vector<RealVec>& r_ij = _r_ij;
...@@ -242,15 +232,14 @@ void ReferenceCCMAAlgorithm::applyConstraints(vector<RealVec>& atomCoordinates, ...@@ -242,15 +232,14 @@ void ReferenceCCMAAlgorithm::applyConstraints(vector<RealVec>& atomCoordinates,
// setup: r_ij for each (i,j) constraint // setup: r_ij for each (i,j) constraint
RealOpenMM tolerance = getTolerance()*2;
for (int ii = 0; ii < _numberOfConstraints; ii++) { for (int ii = 0; ii < _numberOfConstraints; ii++) {
int atomI = _atomIndices[ii].first; int atomI = _atomIndices[ii].first;
int atomJ = _atomIndices[ii].second; int atomJ = _atomIndices[ii].second;
r_ij[ii] = atomCoordinates[atomI] - atomCoordinates[atomJ]; r_ij[ii] = atomCoordinates[atomI] - atomCoordinates[atomJ];
d_ij2[ii] = r_ij[ii].dot(r_ij[ii]); d_ij2[ii] = r_ij[ii].dot(r_ij[ii]);
} }
RealOpenMM lowerTol = 1-2*getTolerance()+getTolerance()*getTolerance(); RealOpenMM lowerTol = 1-2*tolerance+tolerance*tolerance;
RealOpenMM upperTol = 1+2*getTolerance()+getTolerance()*getTolerance(); RealOpenMM upperTol = 1+2*tolerance+tolerance*tolerance;
// main loop // main loop
...@@ -267,7 +256,7 @@ void ReferenceCCMAAlgorithm::applyConstraints(vector<RealVec>& atomCoordinates, ...@@ -267,7 +256,7 @@ void ReferenceCCMAAlgorithm::applyConstraints(vector<RealVec>& atomCoordinates,
if (constrainingVelocities) { if (constrainingVelocities) {
RealOpenMM rrpr = rp_ij.dot(r_ij[ii]); RealOpenMM rrpr = rp_ij.dot(r_ij[ii]);
constraintDelta[ii] = -2*reducedMasses[ii]*rrpr/d_ij2[ii]; constraintDelta[ii] = -2*reducedMasses[ii]*rrpr/d_ij2[ii];
if (fabs(constraintDelta[ii]) <= getTolerance()) if (fabs(constraintDelta[ii]) <= tolerance)
numberConverged++; numberConverged++;
} }
else { else {
......
...@@ -39,8 +39,7 @@ ...@@ -39,8 +39,7 @@
using namespace OpenMM; using namespace OpenMM;
using namespace std; using namespace std;
ReferenceConstraints::ReferenceConstraints(const System& system, RealOpenMM tolerance) : ccma(NULL), settle(NULL) { ReferenceConstraints::ReferenceConstraints(const System& system) : ccma(NULL), settle(NULL) {
this->tolerance = tolerance;
int numParticles = system.getNumParticles(); int numParticles = system.getNumParticles();
vector<RealOpenMM> masses(numParticles); vector<RealOpenMM> masses(numParticles);
for (int i = 0; i < numParticles; ++i) for (int i = 0; i < numParticles; ++i)
...@@ -141,7 +140,7 @@ ReferenceConstraints::ReferenceConstraints(const System& system, RealOpenMM tole ...@@ -141,7 +140,7 @@ ReferenceConstraints::ReferenceConstraints(const System& system, RealOpenMM tole
isSettleAtom[p2] = true; isSettleAtom[p2] = true;
isSettleAtom[p3] = true; isSettleAtom[p3] = true;
} }
settle = new ReferenceSETTLEAlgorithm(atom1, atom2, atom3, distance1, distance2, masses, tolerance); settle = new ReferenceSETTLEAlgorithm(atom1, atom2, atom3, distance1, distance2, masses);
} }
// All other constraints are handled with CCMA. // All other constraints are handled with CCMA.
...@@ -179,7 +178,7 @@ ReferenceConstraints::ReferenceConstraints(const System& system, RealOpenMM tole ...@@ -179,7 +178,7 @@ ReferenceConstraints::ReferenceConstraints(const System& system, RealOpenMM tole
// Create the CCMA object. // Create the CCMA object.
ccma = new ReferenceCCMAAlgorithm(numParticles, numCCMA, ccmaIndices, ccmaDistance, masses, angles, tolerance); ccma = new ReferenceCCMAAlgorithm(numParticles, numCCMA, ccmaIndices, ccmaDistance, masses, angles);
} }
} }
...@@ -190,28 +189,16 @@ ReferenceConstraints::~ReferenceConstraints() { ...@@ -190,28 +189,16 @@ ReferenceConstraints::~ReferenceConstraints() {
delete settle; delete settle;
} }
RealOpenMM ReferenceConstraints::getTolerance() const { void ReferenceConstraints::apply(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) {
return tolerance;
}
void ReferenceConstraints::setTolerance(RealOpenMM tolerance) {
this->tolerance = tolerance;
if (ccma != NULL)
ccma->setTolerance(tolerance);
if (settle != NULL)
settle->setTolerance(tolerance);
}
void ReferenceConstraints::apply(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses) {
if (ccma != NULL) if (ccma != NULL)
ccma->apply(atomCoordinates, atomCoordinatesP, inverseMasses); ccma->apply(atomCoordinates, atomCoordinatesP, inverseMasses, tolerance);
if (settle != NULL) if (settle != NULL)
settle->apply(atomCoordinates, atomCoordinatesP, inverseMasses); settle->apply(atomCoordinates, atomCoordinatesP, inverseMasses, tolerance);
} }
void ReferenceConstraints::applyToVelocities(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& velocities, vector<RealOpenMM>& inverseMasses) { void ReferenceConstraints::applyToVelocities(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& velocities, vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) {
if (ccma != NULL) if (ccma != NULL)
ccma->applyToVelocities(atomCoordinates, velocities, inverseMasses); ccma->applyToVelocities(atomCoordinates, velocities, inverseMasses, tolerance);
if (settle != NULL) if (settle != NULL)
settle->applyToVelocities(atomCoordinates, velocities, inverseMasses); settle->applyToVelocities(atomCoordinates, velocities, inverseMasses, tolerance);
} }
...@@ -97,7 +97,7 @@ ReferenceCustomDynamics::~ReferenceCustomDynamics() { ...@@ -97,7 +97,7 @@ ReferenceCustomDynamics::~ReferenceCustomDynamics() {
void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, vector<RealVec>& atomCoordinates, void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses, vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses,
map<string, RealOpenMM>& globals, vector<vector<RealVec> >& perDof, bool& forcesAreValid){ map<string, RealOpenMM>& globals, vector<vector<RealVec> >& perDof, bool& forcesAreValid, RealOpenMM tolerance){
int numSteps = stepType.size(); int numSteps = stepType.size();
globals.insert(context.getParameters().begin(), context.getParameters().end()); globals.insert(context.getParameters().begin(), context.getParameters().end());
oldPos = atomCoordinates; oldPos = atomCoordinates;
...@@ -252,12 +252,12 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve ...@@ -252,12 +252,12 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
break; break;
} }
case CustomIntegrator::ConstrainPositions: { case CustomIntegrator::ConstrainPositions: {
getReferenceConstraintAlgorithm()->apply(oldPos, atomCoordinates, inverseMasses); getReferenceConstraintAlgorithm()->apply(oldPos, atomCoordinates, inverseMasses, tolerance);
oldPos = atomCoordinates; oldPos = atomCoordinates;
break; break;
} }
case CustomIntegrator::ConstrainVelocities: { case CustomIntegrator::ConstrainVelocities: {
getReferenceConstraintAlgorithm()->applyToVelocities(oldPos, velocities, inverseMasses); getReferenceConstraintAlgorithm()->applyToVelocities(oldPos, velocities, inverseMasses, tolerance);
break; break;
} }
case CustomIntegrator::UpdateContextState: { case CustomIntegrator::UpdateContextState: {
......
...@@ -249,11 +249,12 @@ void ReferenceDynamics::setReferenceConstraintAlgorithm( ReferenceConstraintAlgo ...@@ -249,11 +249,12 @@ void ReferenceDynamics::setReferenceConstraintAlgorithm( ReferenceConstraintAlgo
@param velocities velocities @param velocities velocities
@param forces forces @param forces forces
@param masses atom masses @param masses atom masses
@param tolerance the constraint tolerance
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void ReferenceDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates, void ReferenceDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses ){ vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM tolerance) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
......
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