Skip to content

GitLab

Commit e68df534 authored by peastman's avatar peastman
Browse files

Lots of code cleanup related to constraints in the reference platform

parent 3e36fd7e
Show whitespace changes
Inline Side-by-side
Showing with 366 additions and 1628 deletions
platforms/reference/include/ReferenceCCMAAlgorithm.h
View file @ e68df534
/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -34,7 +34,7 @@
class OPENMM_EXPORT ReferenceCCMAAlgorithm : public ReferenceConstraintAlgorithm {
protected:
protected:
int _maximumNumberOfIterations;
RealOpenMM _tolerance;
......@@ -50,135 +50,73 @@ class OPENMM_EXPORT ReferenceCCMAAlgorithm : public ReferenceConstraintAlgorithm
bool _hasInitializedMasses;
std::vector<std::vector<std::pair<int, RealOpenMM> > > _matrix;
private:
private:
int applyConstraints(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
void applyConstraints(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses, bool constrainingVelocities);
public:
public:
class AngleInfo;
/**---------------------------------------------------------------------------------------
ReferenceCCMAAlgorithm constructor
@param numberOfAtoms number of atoms
@param numberOfConstraints number of constraints
@param atomIndices atom indices for contraints
@param distance distances for constraints
@param masses atom masses
@param angles angle force field terms
@param tolerance constraint tolerance
--------------------------------------------------------------------------------------- */
/**
* Create a ReferenceCCMAAlgorithm object.
*
* @param numberOfAtoms the number of atoms in the system
* @param numberOfConstraints the number of constraints
* @param atomIndices atom indices for contraints
* @param distance distances for constraints
* @param masses atom masses
* @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 );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceCCMAAlgorithm( );
/**---------------------------------------------------------------------------------------
Get number of constraints
@return number of constraints
--------------------------------------------------------------------------------------- */
/**
* Get the number of constraints.
*/
int getNumberOfConstraints( void ) const;
/**---------------------------------------------------------------------------------------
Get maximum number of iterations
@return maximum number of iterations
--------------------------------------------------------------------------------------- */
/**
* Get the maximum number of iterations to perform.
*/
int getMaximumNumberOfIterations( void ) const;
/**---------------------------------------------------------------------------------------
Set maximum number of iterations
@param maximumNumberOfIterations new maximum number of iterations
--------------------------------------------------------------------------------------- */
/**
* Set the maximum number of iterations to perform.
*/
void setMaximumNumberOfIterations( int maximumNumberOfIterations );
/**---------------------------------------------------------------------------------------
Get tolerance
@return tolerance
--------------------------------------------------------------------------------------- */
/**
* Get the constraint tolerance.
*/
RealOpenMM getTolerance( void ) const;
/**---------------------------------------------------------------------------------------
Set tolerance
@param tolerance new tolerance
--------------------------------------------------------------------------------------- */
/**
* Set the constraint tolerance.
*/
void setTolerance( RealOpenMM tolerance );
/**---------------------------------------------------------------------------------------
Apply CCMA algorithm
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomCoordinatesP atom coordinates prime
@param inverseMasses 1/mass
@return SimTKOpenMMCommon::DefaultReturn if converge; else
return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int apply( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses );
/**---------------------------------------------------------------------------------------
Apply constraint algorithm to velocities.
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities atom velocities
@param inverseMasses 1/mass
@return SimTKOpenMMCommon::DefaultReturn if converge; else
return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
/**
* Apply the constraint algorithm.
*
* @param atomCoordinates the original atom coordinates
* @param atomCoordinatesP the new atom coordinates
* @param inverseMasses 1/mass
*/
void apply(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses);
int applyToVelocities(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
/**
* Apply the constraint algorithm to velocities.
*
* @param atomCoordinates the atom coordinates
* @param atomCoordinatesP the velocities to modify
* @param inverseMasses 1/mass
*/
void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses);
/**---------------------------------------------------------------------------------------
Report any violated constraints
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param message report
@return number of violated constraints
--------------------------------------------------------------------------------------- */
int reportCCMA( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::stringstream& message );
};
class ReferenceCCMAAlgorithm::AngleInfo
......
platforms/reference/include/ReferenceConstraint.h deleted 100644 → 0
View file @ 3e36fd7e
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __ReferenceConstraint_H__
#define __ReferenceConstraint_H__
// ---------------------------------------------------------------------------------------
class ReferenceConstraint {
private:
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
ReferenceConstraint( );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceConstraint( );
};
class ReferenceShakeConstraint : public ReferenceConstraint {
private:
// force heavy atom into index 0
int _atomIndices[2];
RealOpenMM _inverseMasses[2];
RealOpenMM _constraintDistance;
public:
/**---------------------------------------------------------------------------------------
Constructor
@param atomIndex1 atom index 1
@param atomIndex2 atom index 2
@param distance distance constraint
@param inverseMass1 inverseMass of atom 1
@param inverseMass2 inverseMass of atom 2
--------------------------------------------------------------------------------------- */
ReferenceShakeConstraint( int atomIndex1, int atomIndex2, RealOpenMM distance,
RealOpenMM inverseMass1, RealOpenMM inverseMass2 );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceShakeConstraint( );
/**---------------------------------------------------------------------------------------
Get constraint distance
@return distance
--------------------------------------------------------------------------------------- */
RealOpenMM getConstraintDistance( void );
/**---------------------------------------------------------------------------------------
Get inverse mass of heavy atom
@return inverse mass
--------------------------------------------------------------------------------------- */
RealOpenMM getHeavyAtomInverseMass( void );
/**---------------------------------------------------------------------------------------
Get inverse mass of light atom
@return inverse mass
--------------------------------------------------------------------------------------- */
RealOpenMM getLightAtomInverseMass( void );
/**---------------------------------------------------------------------------------------
Get index of heavy atom
@return index
--------------------------------------------------------------------------------------- */
int getHeavyAtomIndex( void );
/**---------------------------------------------------------------------------------------
Get index of light atom
@return index
--------------------------------------------------------------------------------------- */
int getLightAtomIndex( void );
};
// ---------------------------------------------------------------------------------------
#endif // __ReferenceShakeConstraint_H__
platforms/reference/include/ReferenceConstraintAlgorithm.h
View file @ e68df534
/* Portions copyright (c) 2009 Stanford University and Simbios.
/* Portions copyright (c) 2009-2013 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -36,54 +36,34 @@ public:
virtual ~ReferenceConstraintAlgorithm() {};
/**---------------------------------------------------------------------------------------
Get the constraint tolerance
--------------------------------------------------------------------------------------- */
/**
* Get the constraint tolerance.
*/
virtual RealOpenMM getTolerance() const = 0;
/**---------------------------------------------------------------------------------------
Set the constraint tolerance
--------------------------------------------------------------------------------------- */
/**
* Set the constraint tolerance.
*/
virtual void setTolerance(RealOpenMM tolerance) = 0;
/**---------------------------------------------------------------------------------------
Apply constraint algorithm
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomCoordinatesP atom coordinates prime
@param inverseMasses 1/mass
@return SimTKOpenMMCommon::DefaultReturn if converge; else
return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
virtual int apply(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
/**
* Apply the constraint algorithm.
*
* @param atomCoordinates the original atom coordinates
* @param atomCoordinatesP the new atom coordinates
* @param inverseMasses 1/mass
*/
virtual void apply(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses) = 0;
/**---------------------------------------------------------------------------------------
Apply constraint algorithm to velocities.
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities atom velocities
@param inverseMasses 1/mass
@return SimTKOpenMMCommon::DefaultReturn if converge; else
return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
virtual int applyToVelocities(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
/**
* Apply the constraint algorithm to velocities.
*
* @param atomCoordinates the atom coordinates
* @param atomCoordinatesP the velocities to modify
* @param inverseMasses 1/mass
*/
virtual void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses) = 0;
};
......
platforms/reference/include/ReferenceConstraints.h
View file @ e68df534
......@@ -61,22 +61,20 @@ public:
/**
* Apply the constraint algorithm.
*
* @param numberOfAtoms number of atoms
* @param atomCoordinates the original atom coordinates
* @param atomCoordinatesP the new atom coordinates
* @param inverseMasses 1/mass
*/
int apply(int numberOfAtoms, 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);
/**
* Apply the constraint algorithm to velocities.
*
* @param numberOfAtoms number of atoms
* @param atomCoordinates the atom coordinates
* @param atomCoordinatesP the velocities to modify
* @param inverseMasses 1/mass
*/
int applyToVelocities(int numberOfAtoms, 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);
private:
RealOpenMM tolerance;
ReferenceCCMAAlgorithm* ccma;
......
platforms/reference/include/ReferenceSETTLEAlgorithm.h
View file @ e68df534
......@@ -58,22 +58,20 @@ public:
/**
* Apply the constraint algorithm.
*
* @param numberOfAtoms number of atoms
* @param atomCoordinates the original atom coordinates
* @param atomCoordinatesP the new atom coordinates
* @param inverseMasses 1/mass
*/
int apply(int numberOfAtoms, 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);
/**
* Apply the constraint algorithm to velocities.
*
* @param numberOfAtoms number of atoms
* @param atomCoordinates the atom coordinates
* @param atomCoordinatesP the velocities to modify
* @param inverseMasses 1/mass
*/
int applyToVelocities(int numberOfAtoms, 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);
private:
RealOpenMM tolerance;
std::vector<int> atom1;
......
platforms/reference/include/ReferenceShakeAlgorithm.h deleted 100644 → 0
View file @ 3e36fd7e
/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __ReferenceShakeAlgorithm_H__
#define __ReferenceShakeAlgorithm_H__
#include "ReferenceConstraintAlgorithm.h"
// ---------------------------------------------------------------------------------------
class ReferenceShakeAlgorithm : public ReferenceConstraintAlgorithm {
protected:
int _maximumNumberOfIterations;
RealOpenMM _tolerance;
int _numberOfConstraints;
int** _atomIndices;
RealOpenMM* _distance;
RealOpenMM** _r_ij;
RealOpenMM* _d_ij2;
RealOpenMM* _distanceTolerance;
RealOpenMM* _reducedMasses;
bool _hasInitializedMasses;
public:
/**---------------------------------------------------------------------------------------
ReferenceShakeAlgorithm constructor
@param numberOfConstraints number of constraints
@param atomIndices atom indices for contraints
@param distance distances for constraints
@param tolerance constraint tolerance
--------------------------------------------------------------------------------------- */
ReferenceShakeAlgorithm( int numberOfConstraints, int** atomIndices, RealOpenMM* distance, RealOpenMM tolerance );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceShakeAlgorithm( );
/**---------------------------------------------------------------------------------------
Get number of constraints
@return number of constraints
--------------------------------------------------------------------------------------- */
int getNumberOfConstraints( void ) const;
/**---------------------------------------------------------------------------------------
Get maximum number of iterations
@return maximum number of iterations
--------------------------------------------------------------------------------------- */
int getMaximumNumberOfIterations( void ) const;
/**---------------------------------------------------------------------------------------
Set maximum number of iterations
@param maximumNumberOfIterations new maximum number of iterations
--------------------------------------------------------------------------------------- */
void setMaximumNumberOfIterations( int maximumNumberOfIterations );
/**---------------------------------------------------------------------------------------
Get tolerance
@return tolerance
--------------------------------------------------------------------------------------- */
RealOpenMM getTolerance( void ) const;
/**---------------------------------------------------------------------------------------
Set tolerance
@param tolerance new tolerance
--------------------------------------------------------------------------------------- */
void setTolerance( RealOpenMM tolerance );
/**---------------------------------------------------------------------------------------
Apply Shake algorithm
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomCoordinatesP atom coordinates prime
@param inverseMasses 1/mass
@return SimTKOpenMMCommon::DefaultReturn if converge; else
return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int apply( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& atomCoordinatesP, std::vector<RealOpenMM>& inverseMasses );
/**---------------------------------------------------------------------------------------
Apply constraint algorithm to velocities.
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities atom velocities
@param inverseMasses 1/mass
@return SimTKOpenMMCommon::DefaultReturn if converge; else
return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int applyToVelocities(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses);
/**---------------------------------------------------------------------------------------
Report any violated constraints
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param message report
@return number of violated constraints
--------------------------------------------------------------------------------------- */
int reportShake( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, std::stringstream& message );
};
// ---------------------------------------------------------------------------------------
#endif // __ReferenceShakeAlgorithm_H__
platforms/reference/src/ReferenceKernels.cpp
View file @ e68df534
......@@ -160,7 +160,7 @@ static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>&
for (int i = 0; i < numParticles; i++)
inverseMasses[i] = (masses[i] == 0 ? 0 : 1/masses[i]);
constraints->setTolerance(1e-4);
constraints->applyToVelocities(numParticles, posData, shiftedVel, inverseMasses);
constraints->applyToVelocities(posData, shiftedVel, inverseMasses);
}
// Compute the kinetic energy.
......@@ -324,7 +324,7 @@ void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
constraints = new ReferenceConstraints(context.getSystem(), (RealOpenMM) tol);
vector<RealVec>& positions = extractPositions(context);
constraints->setTolerance(tol);
constraints->apply(data.numParticles, positions, positions, inverseMasses);
constraints->apply(positions, positions, inverseMasses);
ReferenceVirtualSites::computePositions(context.getSystem(), positions);
}
......@@ -334,7 +334,7 @@ void ReferenceApplyConstraintsKernel::applyToVelocities(ContextImpl& context, do
vector<RealVec>& positions = extractPositions(context);
vector<RealVec>& velocities = extractVelocities(context);
constraints->setTolerance(tol);
constraints->applyToVelocities(data.numParticles, positions, velocities, inverseMasses);
constraints->applyToVelocities(positions, velocities, inverseMasses);
}
void ReferenceVirtualSitesKernel::initialize(const System& system) {
......
platforms/reference/src/SimTKReference/ReferenceBrownianDynamics.cpp
View file @ e68df534
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* Contributors: Peter Eastman, Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -161,8 +161,8 @@ void ReferenceBrownianDynamics::update(const OpenMM::System& system, vector<Real
}
}
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if( referenceConstraintAlgorithm )
referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime, inverseMasses );
if (referenceConstraintAlgorithm)
referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses);
// Update the positions and velocities.
......
platforms/reference/src/SimTKReference/ReferenceCCMAAlgorithm.cpp
View file @ e68df534
/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* Contributors: Peter Eastman, Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -31,6 +31,7 @@
#include "ReferenceCCMAAlgorithm.h"
#include "ReferenceDynamics.h"
#include "quern.h"
#include "openmm/OpenMMException.h"
#include "openmm/Vec3.h"
#include <map>
......@@ -38,45 +39,17 @@ using std::map;
using std::pair;
using std::vector;
using std::set;
using OpenMM::OpenMMException;
using OpenMM::Vec3;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
ReferenceCCMAAlgorithm constructor
@param numberOfAtoms number of atoms
@param numberOfConstraints number of constraints
@param atomIndices atom indices for contraints
@param distance distances for constraints
@param masses atom masses
@param angles angle force field terms
@param tolerance constraint tolerance
--------------------------------------------------------------------------------------- */
ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm( int numberOfAtoms,
ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm(int numberOfAtoms,
int numberOfConstraints,
const vector<pair<int, int> >& atomIndices,
const vector<RealOpenMM>& distance,
vector<RealOpenMM>& masses,
vector<AngleInfo>& angles,
RealOpenMM tolerance){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM oneM = -1.0;
static const int threeI = 3;
// ---------------------------------------------------------------------------------------
RealOpenMM tolerance) {
_numberOfConstraints = numberOfConstraints;
_atomIndices = atomIndices;
_distance = distance;
......@@ -89,9 +62,9 @@ ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm( int numberOfAtoms,
if (_numberOfConstraints > 0) {
_r_ij.resize(numberOfConstraints);
_d_ij2 = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "dij_2" );
_distanceTolerance = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "distanceTolerance" );
_reducedMasses = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "reducedMasses" );
_d_ij2 = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray(numberOfConstraints, NULL, 1, 0.0, "dij_2");
_distanceTolerance = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray(numberOfConstraints, NULL, 1, 0.0, "distanceTolerance");
_reducedMasses = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray(numberOfConstraints, NULL, 1, 0.0, "reducedMasses");
}
if (numberOfConstraints > 0)
{
......@@ -112,8 +85,8 @@ ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm( int numberOfAtoms,
int atomj1 = _atomIndices[j].second;
int atomk0 = _atomIndices[k].first;
int atomk1 = _atomIndices[k].second;
RealOpenMM invMass0 = one/masses[atomj0];
RealOpenMM invMass1 = one/masses[atomj1];
RealOpenMM invMass0 = 1/masses[atomj0];
RealOpenMM invMass1 = 1/masses[atomj1];
int atoma, atomb, atomc;
if (atomj0 == atomk0) {
atoma = atomj1;
......@@ -208,178 +181,48 @@ ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm( int numberOfAtoms,
}
}
/**---------------------------------------------------------------------------------------
ReferenceCCMAAlgorithm destructor
--------------------------------------------------------------------------------------- */
ReferenceCCMAAlgorithm::~ReferenceCCMAAlgorithm( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCCMAAlgorithm::~ReferenceCCMAAlgorithm";
// ---------------------------------------------------------------------------------------
ReferenceCCMAAlgorithm::~ReferenceCCMAAlgorithm() {
if (_numberOfConstraints > 0) {
SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _d_ij2, "d_ij2" );
SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _distanceTolerance, "distanceTolerance" );
SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _reducedMasses, "reducedMasses" );
SimTKOpenMMUtilities::freeOneDRealOpenMMArray(_d_ij2, "d_ij2");
SimTKOpenMMUtilities::freeOneDRealOpenMMArray(_distanceTolerance, "distanceTolerance");
SimTKOpenMMUtilities::freeOneDRealOpenMMArray(_reducedMasses, "reducedMasses");
}
}
/**---------------------------------------------------------------------------------------
Get number of constraints
@return number of constraints
--------------------------------------------------------------------------------------- */
int ReferenceCCMAAlgorithm::getNumberOfConstraints( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCCMAAlgorithm::getNumberOfConstraints";
// ---------------------------------------------------------------------------------------
int ReferenceCCMAAlgorithm::getNumberOfConstraints() const {
return _numberOfConstraints;
}
/**---------------------------------------------------------------------------------------
Get maximum number of iterations
@return maximum number of iterations
--------------------------------------------------------------------------------------- */
int ReferenceCCMAAlgorithm::getMaximumNumberOfIterations( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCCMAAlgorithm::getMaximumNumberOfIterations";
// ---------------------------------------------------------------------------------------
int ReferenceCCMAAlgorithm::getMaximumNumberOfIterations() const {
return _maximumNumberOfIterations;
}
/**---------------------------------------------------------------------------------------
Set maximum number of iterations
@param maximumNumberOfIterations new maximum number of iterations
--------------------------------------------------------------------------------------- */
void ReferenceCCMAAlgorithm::setMaximumNumberOfIterations( int maximumNumberOfIterations ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCCMAAlgorithm::setMaximumNumberOfIterations";
// ---------------------------------------------------------------------------------------
void ReferenceCCMAAlgorithm::setMaximumNumberOfIterations(int maximumNumberOfIterations) {
_maximumNumberOfIterations = maximumNumberOfIterations;
}
/**---------------------------------------------------------------------------------------
Get tolerance
@return tolerance
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceCCMAAlgorithm::getTolerance( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCCMAAlgorithm::getTolerance";
// ---------------------------------------------------------------------------------------
RealOpenMM ReferenceCCMAAlgorithm::getTolerance() const {
return _tolerance;
}
/**---------------------------------------------------------------------------------------
Set tolerance
@param tolerance new tolerance
--------------------------------------------------------------------------------------- */
void ReferenceCCMAAlgorithm::setTolerance( RealOpenMM tolerance ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCCMAAlgorithm::setTolerance";
// ---------------------------------------------------------------------------------------
_tolerance = tolerance;;
void ReferenceCCMAAlgorithm::setTolerance(RealOpenMM tolerance) {
_tolerance = tolerance;
}
/**---------------------------------------------------------------------------------------
Apply CCMA algorithm
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomCoordinatesP atom coordinates prime
@param inverseMasses 1/mass
@return SimTKOpenMMCommon::DefaultReturn if converge; else
return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int ReferenceCCMAAlgorithm::apply( int numberOfAtoms, vector<RealVec>& atomCoordinates,
void ReferenceCCMAAlgorithm::apply(vector<RealVec>& atomCoordinates,
vector<RealVec>& atomCoordinatesP,
vector<RealOpenMM>& inverseMasses ){
return applyConstraints(numberOfAtoms, atomCoordinates, atomCoordinatesP, inverseMasses, false);
vector<RealOpenMM>& inverseMasses) {
applyConstraints(atomCoordinates, atomCoordinatesP, inverseMasses, false);
}
/**---------------------------------------------------------------------------------------
Apply constraint algorithm to velocities.
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities atom velocities
@param inverseMasses 1/mass
@return SimTKOpenMMCommon::DefaultReturn if converge; else
return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int ReferenceCCMAAlgorithm::applyToVelocities(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
void ReferenceCCMAAlgorithm::applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses) {
return applyConstraints(numberOfAtoms, atomCoordinates, velocities, inverseMasses, true);
applyConstraints(atomCoordinates, velocities, inverseMasses, true);
}
int ReferenceCCMAAlgorithm::applyConstraints(int numberOfAtoms, vector<RealVec>& atomCoordinates,
void ReferenceCCMAAlgorithm::applyConstraints(vector<RealVec>& atomCoordinates,
vector<RealVec>& atomCoordinatesP,
vector<RealOpenMM>& inverseMasses, bool constrainingVelocities){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceCCMAAlgorithm::apply";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM half = 0.5;
static const RealOpenMM epsilon6 = (RealOpenMM) 1.0e-06;
// ---------------------------------------------------------------------------------------
vector<RealOpenMM>& inverseMasses, bool constrainingVelocities) {
// temp arrays
vector<RealVec>& r_ij = _r_ij;
......@@ -388,27 +231,26 @@ int ReferenceCCMAAlgorithm::applyConstraints(int numberOfAtoms, vector<RealVec>&
// calculate reduced masses on 1st pass
if( !_hasInitializedMasses ){
if (!_hasInitializedMasses) {
_hasInitializedMasses = true;
for( int ii = 0; ii < _numberOfConstraints; ii++ ){
for (int ii = 0; ii < _numberOfConstraints; ii++) {
int atomI = _atomIndices[ii].first;
int atomJ = _atomIndices[ii].second;
reducedMasses[ii] = half/( inverseMasses[atomI] + inverseMasses[atomJ] );
reducedMasses[ii] = 0.5/(inverseMasses[atomI] + inverseMasses[atomJ]);
}
}
// setup: r_ij for each (i,j) constraint
RealOpenMM tolerance = getTolerance();
tolerance *= two;
for( int ii = 0; ii < _numberOfConstraints; ii++ ){
RealOpenMM tolerance = getTolerance()*2;
for (int ii = 0; ii < _numberOfConstraints; ii++) {
int atomI = _atomIndices[ii].first;
int atomJ = _atomIndices[ii].second;
r_ij[ii] = atomCoordinates[atomI] - atomCoordinates[atomJ];
d_ij2[ii] = r_ij[ii].dot(r_ij[ii]);
}
RealOpenMM lowerTol = one-two*getTolerance()+getTolerance()*getTolerance();
RealOpenMM upperTol = one+two*getTolerance()+getTolerance()*getTolerance();
RealOpenMM lowerTol = 1-2*getTolerance()+getTolerance()*getTolerance();
RealOpenMM upperTol = 1+2*getTolerance()+getTolerance()*getTolerance();
// main loop
......@@ -418,11 +260,9 @@ int ReferenceCCMAAlgorithm::applyConstraints(int numberOfAtoms, vector<RealVec>&
vector<RealOpenMM> tempDelta(_numberOfConstraints);
while (iterations < getMaximumNumberOfIterations()) {
numberConverged = 0;
for( int ii = 0; ii < _numberOfConstraints; ii++ ){
for (int ii = 0; ii < _numberOfConstraints; ii++) {
int atomI = _atomIndices[ii].first;
int atomJ = _atomIndices[ii].second;
RealVec rp_ij = atomCoordinatesP[atomI] - atomCoordinatesP[atomJ];
if (constrainingVelocities) {
RealOpenMM rrpr = rp_ij.dot(r_ij[ii]);
......@@ -434,20 +274,14 @@ int ReferenceCCMAAlgorithm::applyConstraints(int numberOfAtoms, vector<RealVec>&
RealOpenMM rp2 = rp_ij.dot(rp_ij);
RealOpenMM dist2 = _distance[ii]*_distance[ii];
RealOpenMM diff = dist2 - rp2;
constraintDelta[ii] = zero;
RealOpenMM rrpr = DOT3( rp_ij, r_ij[ii] );
if( rrpr < d_ij2[ii]*epsilon6 ){
std::stringstream message;
message << iterations <<" "<<atomI<<" "<<atomJ<< " Error: sign of rrpr < 0?\n";
SimTKOpenMMLog::printMessage( message );
} else {
constraintDelta[ii] = 0;
RealOpenMM rrpr = DOT3(rp_ij, r_ij[ii]);
constraintDelta[ii] = reducedMasses[ii]*diff/rrpr;
}
if (rp2 >= lowerTol*dist2 && rp2 <= upperTol*dist2)
numberConverged++;
}
}
if( numberConverged == _numberOfConstraints )
if (numberConverged == _numberOfConstraints)
break;
iterations++;
......@@ -462,8 +296,7 @@ int ReferenceCCMAAlgorithm::applyConstraints(int numberOfAtoms, vector<RealVec>&
}
constraintDelta = tempDelta;
}
for( int ii = 0; ii < _numberOfConstraints; ii++ ){
for (int ii = 0; ii < _numberOfConstraints; ii++) {
int atomI = _atomIndices[ii].first;
int atomJ = _atomIndices[ii].second;
RealVec dr = r_ij[ii]*constraintDelta[ii];
......@@ -471,68 +304,4 @@ int ReferenceCCMAAlgorithm::applyConstraints(int numberOfAtoms, vector<RealVec>&
atomCoordinatesP[atomJ] -= dr*inverseMasses[atomJ];
}
}
return (numberConverged == _numberOfConstraints ? SimTKOpenMMCommon::DefaultReturn : SimTKOpenMMCommon::ErrorReturn);
}
/**---------------------------------------------------------------------------------------
Report any violated constriants
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param message report
@return number of violated constraints
--------------------------------------------------------------------------------------- */
int ReferenceCCMAAlgorithm::reportCCMA( int numberOfAtoms, vector<RealVec>& atomCoordinates,
std::stringstream& message ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceCCMAAlgorithm::reportCCMA";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM oneM = -1.0;
static const RealOpenMM half = 0.5;
// ---------------------------------------------------------------------------------------
int numberOfConstraints = getNumberOfConstraints();
// loop over constraints calculating distance and comparing to
// expected distance -- report any contraints that are violated
int numberConverged = 0;
RealOpenMM tolerance = getTolerance();
for( int ii = 0; ii < _numberOfConstraints; ii++ ){
int atomI = _atomIndices[ii].first;
int atomJ = _atomIndices[ii].second;
RealOpenMM rp2 = zero;
for( int jj = 0; jj < 3; jj++ ){
rp2 += (atomCoordinates[atomI][jj] - atomCoordinates[atomJ][jj])*(atomCoordinates[atomI][jj] - atomCoordinates[atomJ][jj]);
}
RealOpenMM diff = FABS( rp2 - (_distance[ii]*_distance[ii]) );
if( diff > tolerance ){
message << ii << " constraint violated: " << atomI << " " << atomJ << "] d=" << SQRT( rp2 ) << " " << rp2 << " d0=" << _distance[ii];
message << " diff=" << diff;
message << " [" << atomCoordinates[atomI][0] << " " << atomCoordinates[atomI][1] << " " << atomCoordinates[atomI][2] << "] ";
message << " [" << atomCoordinates[atomJ][0] << " " << atomCoordinates[atomJ][1] << " " << atomCoordinates[atomJ][2] << "] ";
message << "\n";
} else {
numberConverged++;
}
}
return (numberOfConstraints-numberConverged);
}
platforms/reference/src/SimTKReference/ReferenceConstraint.cpp deleted 100644 → 0
View file @ 3e36fd7e
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* 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 <string.h>
#include <sstream>
#include "SimTKOpenMMCommon.h"
#include "SimTKOpenMMLog.h"
#include "SimTKOpenMMUtilities.h"
#include "ReferenceConstraint.h"
#include "ReferenceDynamics.h"
/**---------------------------------------------------------------------------------------
ReferenceConstraint constructor
--------------------------------------------------------------------------------------- */
ReferenceConstraint::ReferenceConstraint( ){
// ---------------------------------------------------------------------------------------
//static const char* methodName = "\nReferenceConstraint::ReferenceConstraint";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
ReferenceConstraint destructor
--------------------------------------------------------------------------------------- */
ReferenceConstraint::~ReferenceConstraint( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceConstraint::~ReferenceConstraint";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
ReferenceShakeConstraint constructor
@param numberOfAtoms number of atoms
@param deltaT delta t for dynamics
@param tau viscosity
@param temperature temperature
--------------------------------------------------------------------------------------- */
ReferenceShakeConstraint::ReferenceShakeConstraint( int atomIndex1, int atomIndex2,
RealOpenMM constraintDistance,
RealOpenMM inverseMass1,
RealOpenMM inverseMass2 ) : ReferenceConstraint( ) {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceShakeConstraint::ReferenceShakeConstraint";
// ---------------------------------------------------------------------------------------
// insure heavy atom is in 0 slot
if( inverseMass1 > inverseMass2 ){
int temp = atomIndex1;
atomIndex1 = atomIndex2;
atomIndex2 = temp;
RealOpenMM tempR = inverseMass1;
inverseMass1 = inverseMass2;
inverseMass2 = tempR;
}
_atomIndices[0] = atomIndex1;
_atomIndices[1] = atomIndex2;
_inverseMasses[0] = inverseMass1;
_inverseMasses[1] = inverseMass2;
_constraintDistance = constraintDistance;
}
/**---------------------------------------------------------------------------------------
ReferenceShakeConstraint destructor
--------------------------------------------------------------------------------------- */
ReferenceShakeConstraint::~ReferenceShakeConstraint( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeConstraint::~ReferenceShakeConstraint";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
Get constraint distance
@return constraintDistance
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceShakeConstraint::getConstraintDistance( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeConstraint::getConstraintDistance";
// ---------------------------------------------------------------------------------------
return _constraintDistance;
}
/**---------------------------------------------------------------------------------------
Get inverse mass of heavy atom
@return inverse mass
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceShakeConstraint::getHeavyAtomInverseMass( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeConstraint::getHeavyAtomInverseMass";
// ---------------------------------------------------------------------------------------
return _inverseMasses[0];
}
/**---------------------------------------------------------------------------------------
Get inverse mass of light atom
@return inverse mass
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceShakeConstraint::getLightAtomInverseMass( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeConstraint::getLightAtomInverseMass";
// ---------------------------------------------------------------------------------------
return _inverseMasses[1];
}
/**---------------------------------------------------------------------------------------
Get index of heavy atom
@return index
--------------------------------------------------------------------------------------- */
int ReferenceShakeConstraint::getHeavyAtomIndex( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeConstraint::getHeavyAtomIndex";
// ---------------------------------------------------------------------------------------
return _atomIndices[0];
}
/**---------------------------------------------------------------------------------------
Get index of light atom
@return index
--------------------------------------------------------------------------------------- */
int ReferenceShakeConstraint::getLightAtomIndex( void ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeConstraint::getLightAtomIndex";
// ---------------------------------------------------------------------------------------
return _atomIndices[1];
}
platforms/reference/src/SimTKReference/ReferenceConstraints.cpp
View file @ e68df534
......@@ -202,18 +202,16 @@ void ReferenceConstraints::setTolerance(RealOpenMM tolerance) {
settle->setTolerance(tolerance);
}
int ReferenceConstraints::apply(int numberOfAtoms, vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses) {
void ReferenceConstraints::apply(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses) {
if (ccma != NULL)
ccma->apply(numberOfAtoms, atomCoordinates, atomCoordinatesP, inverseMasses);
ccma->apply(atomCoordinates, atomCoordinatesP, inverseMasses);
if (settle != NULL)
settle->apply(numberOfAtoms, atomCoordinates, atomCoordinatesP, inverseMasses);
return SimTKOpenMMCommon::DefaultReturn;
settle->apply(atomCoordinates, atomCoordinatesP, inverseMasses);
}
int ReferenceConstraints::applyToVelocities(int numberOfAtoms, 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) {
if (ccma != NULL)
ccma->applyToVelocities(numberOfAtoms, atomCoordinates, velocities, inverseMasses);
ccma->applyToVelocities(atomCoordinates, velocities, inverseMasses);
if (settle != NULL)
settle->applyToVelocities(numberOfAtoms, atomCoordinates, velocities, inverseMasses);
return SimTKOpenMMCommon::DefaultReturn;
settle->applyToVelocities(atomCoordinates, velocities, inverseMasses);
}
platforms/reference/src/SimTKReference/ReferenceCustomDynamics.cpp
View file @ e68df534
/* Portions copyright (c) 2011 Stanford University and Simbios.
/* Portions copyright (c) 2011-2013 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -252,12 +252,12 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
break;
}
case CustomIntegrator::ConstrainPositions: {
getReferenceConstraintAlgorithm()->apply(numberOfAtoms, oldPos, atomCoordinates, inverseMasses);
getReferenceConstraintAlgorithm()->apply(oldPos, atomCoordinates, inverseMasses);
oldPos = atomCoordinates;
break;
}
case CustomIntegrator::ConstrainVelocities: {
getReferenceConstraintAlgorithm()->applyToVelocities(numberOfAtoms, oldPos, velocities, inverseMasses);
getReferenceConstraintAlgorithm()->applyToVelocities(oldPos, velocities, inverseMasses);
break;
}
case CustomIntegrator::UpdateContextState: {
......
platforms/reference/src/SimTKReference/ReferenceSETTLEAlgorithm.cpp
View file @ e68df534
......@@ -47,7 +47,7 @@ void ReferenceSETTLEAlgorithm::setTolerance(RealOpenMM tolerance) {
this->tolerance = tolerance;
}
int ReferenceSETTLEAlgorithm::apply(int numberOfAtoms, vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses) {
void ReferenceSETTLEAlgorithm::apply(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses) {
for (int index = 0; index < (int) atom1.size(); ++index) {
RealVec apos0 = atomCoordinates[atom1[index]];
RealVec xp0 = atomCoordinatesP[atom1[index]]-apos0;
......@@ -188,10 +188,9 @@ int ReferenceSETTLEAlgorithm::apply(int numberOfAtoms, vector<OpenMM::RealVec>&
atomCoordinatesP[atom2[index]] = xp1+apos1;
atomCoordinatesP[atom3[index]] = xp2+apos2;
}
return SimTKOpenMMCommon::DefaultReturn;
}
int ReferenceSETTLEAlgorithm::applyToVelocities(int numberOfAtoms, vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& velocities, vector<RealOpenMM>& inverseMasses) {
void ReferenceSETTLEAlgorithm::applyToVelocities(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& velocities, vector<RealOpenMM>& inverseMasses) {
for (int index = 0; index < (int) atom1.size(); ++index) {
RealVec apos0 = atomCoordinates[atom1[index]];
RealVec apos1 = atomCoordinates[atom2[index]];
......@@ -238,5 +237,4 @@ int ReferenceSETTLEAlgorithm::applyToVelocities(int numberOfAtoms, vector<OpenMM
velocities[atom2[index]] = v1;
velocities[atom3[index]] = v2;
}
return SimTKOpenMMCommon::DefaultReturn;
}
platforms/reference/src/SimTKReference/ReferenceShakeAlgorithm.cpp deleted 100644 → 0
View file @ 3e36fd7e
/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
* Contributors: Pande Group
*
* 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 <string.h>
#include <sstream>
#include "SimTKOpenMMCommon.h"
#include "SimTKOpenMMUtilities.h"
#include "SimTKOpenMMLog.h"
#include "ReferenceShakeAlgorithm.h"
#include "ReferenceDynamics.h"
#include "openmm/OpenMMException.h"
using std::vector;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
ReferenceShakeAlgorithm constructor
@param numberOfConstraints number of constraints
@param atomIndices block of atom indices
@param shakeParameters Shake parameters
@param tolerance constraint tolerance
--------------------------------------------------------------------------------------- */
ReferenceShakeAlgorithm::ReferenceShakeAlgorithm( int numberOfConstraints,
int** atomIndices,
RealOpenMM* distance,
RealOpenMM tolerance){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeAlgorithm::ReferenceShakeAlgorithm";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM oneM = -1.0;
static const int threeI = 3;
// ---------------------------------------------------------------------------------------
_numberOfConstraints = numberOfConstraints;
_atomIndices = atomIndices;
_distance = distance;
_maximumNumberOfIterations = 150;
_tolerance = tolerance;
_hasInitializedMasses = false;
// work arrays
if (_numberOfConstraints > 0) {
_r_ij = SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray( numberOfConstraints, threeI, NULL,
1, zero, "r_ij" );
_d_ij2 = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "dij_2" );
_distanceTolerance = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "distanceTolerance" );
_reducedMasses = SimTKOpenMMUtilities::allocateOneDRealOpenMMArray( numberOfConstraints, NULL, 1, zero, "reducedMasses" );
}
}
/**---------------------------------------------------------------------------------------
ReferenceShakeAlgorithm destructor
--------------------------------------------------------------------------------------- */
ReferenceShakeAlgorithm::~ReferenceShakeAlgorithm( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeAlgorithm::~ReferenceShakeAlgorithm";
// ---------------------------------------------------------------------------------------
if (_numberOfConstraints > 0) {
SimTKOpenMMUtilities::freeTwoDRealOpenMMArray( _r_ij, "r_ij" );
SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _d_ij2, "d_ij2" );
SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _distanceTolerance, "distanceTolerance" );
SimTKOpenMMUtilities::freeOneDRealOpenMMArray( _reducedMasses, "reducedMasses" );
}
}
/**---------------------------------------------------------------------------------------
Get number of constraints
@return number of constraints
--------------------------------------------------------------------------------------- */
int ReferenceShakeAlgorithm::getNumberOfConstraints( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeAlgorithm::getNumberOfConstraints";
// ---------------------------------------------------------------------------------------
return _numberOfConstraints;
}
/**---------------------------------------------------------------------------------------
Get maximum number of iterations
@return maximum number of iterations
--------------------------------------------------------------------------------------- */
int ReferenceShakeAlgorithm::getMaximumNumberOfIterations( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeAlgorithm::getMaximumNumberOfIterations";
// ---------------------------------------------------------------------------------------
return _maximumNumberOfIterations;
}
/**---------------------------------------------------------------------------------------
Set maximum number of iterations
@param maximumNumberOfIterations new maximum number of iterations
--------------------------------------------------------------------------------------- */
void ReferenceShakeAlgorithm::setMaximumNumberOfIterations( int maximumNumberOfIterations ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeAlgorithm::setMaximumNumberOfIterations";
// ---------------------------------------------------------------------------------------
_maximumNumberOfIterations = maximumNumberOfIterations;
}
/**---------------------------------------------------------------------------------------
Get tolerance
@return tolerance
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceShakeAlgorithm::getTolerance( void ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeAlgorithm::getTolerance";
// ---------------------------------------------------------------------------------------
return _tolerance;
}
/**---------------------------------------------------------------------------------------
Set tolerance
@param tolerance new tolerance
--------------------------------------------------------------------------------------- */
void ReferenceShakeAlgorithm::setTolerance( RealOpenMM tolerance ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceShakeAlgorithm::setTolerance";
// ---------------------------------------------------------------------------------------
_tolerance = tolerance;;
}
/**---------------------------------------------------------------------------------------
Apply Shake algorithm
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomCoordinatesP atom coordinates prime
@param inverseMasses 1/mass
@return SimTKOpenMMCommon::DefaultReturn if converge; else
return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int ReferenceShakeAlgorithm::apply( int numberOfAtoms, vector<RealVec>& atomCoordinates,
vector<RealVec>& atomCoordinatesP,
vector<RealOpenMM>& inverseMasses ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceShakeAlgorithm::apply";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM oneM = -1.0;
static const RealOpenMM half = 0.5;
static const RealOpenMM epsilon6 = (RealOpenMM) 1.0e-06;
// ---------------------------------------------------------------------------------------
int numberOfConstraints = getNumberOfConstraints();
// temp arrays
RealOpenMM** r_ij = _r_ij;
RealOpenMM* d_ij2 = _d_ij2;
RealOpenMM* distanceTolerance = _distanceTolerance;
RealOpenMM* reducedMasses = _reducedMasses;
// calculate reduced masses on 1st pass
if( !_hasInitializedMasses ){
_hasInitializedMasses = true;
for( int ii = 0; ii < _numberOfConstraints; ii++ ){
int atomI = _atomIndices[ii][0];
int atomJ = _atomIndices[ii][1];
reducedMasses[ii] = half/( inverseMasses[atomI] + inverseMasses[atomJ] );
}
}
// setup: r_ij for each (i,j) constraint
RealOpenMM tolerance = getTolerance();
tolerance *= two;
for( int ii = 0; ii < _numberOfConstraints; ii++ ){
int atomI = _atomIndices[ii][0];
int atomJ = _atomIndices[ii][1];
for( int jj = 0; jj < 3; jj++ ){
r_ij[ii][jj] = atomCoordinates[atomI][jj] - atomCoordinates[atomJ][jj];
}
d_ij2[ii] = DOT3( r_ij[ii], r_ij[ii] );
distanceTolerance[ii] = d_ij2[ii]*tolerance;
if( distanceTolerance[ii] > zero ){
distanceTolerance[ii] = one/distanceTolerance[ii];
}
}
// main loop
int done = 0;
int iterations = 0;
int numberConverged = 0;
while( !done && iterations++ < getMaximumNumberOfIterations() ){
numberConverged = 0;
for( int ii = 0; ii < _numberOfConstraints; ii++ ){
int atomI = _atomIndices[ii][0];
int atomJ = _atomIndices[ii][1];
RealOpenMM rp_ij[3];
for( int jj = 0; jj < 3; jj++ ){
rp_ij[jj] = atomCoordinatesP[atomI][jj] - atomCoordinatesP[atomJ][jj];
}
RealOpenMM rp2 = DOT3( rp_ij, rp_ij );
RealOpenMM dist2= _distance[ii]*_distance[ii];
RealOpenMM diff = dist2 - rp2;
int iconv = (int) ( FABS( diff )*distanceTolerance[ii] );
if( iconv ){
RealOpenMM rrpr = DOT3( rp_ij, r_ij[ii] );
RealOpenMM acor;
if( rrpr < d_ij2[ii]*epsilon6 ){
std::stringstream message;
message << iterations << " Error: sign of rrpr < 0?\n";
SimTKOpenMMLog::printMessage( message );
} else {
acor = reducedMasses[ii]*diff/rrpr;
for( int jj = 0; jj < 3; jj++ ){
RealOpenMM dr = acor*r_ij[ii][jj];
atomCoordinatesP[atomI][jj] += inverseMasses[atomI]*dr;
atomCoordinatesP[atomJ][jj] -= inverseMasses[atomJ]*dr;
}
}
} else {
numberConverged++;
}
}
if( numberConverged == _numberOfConstraints ){
done = true;
}
}
return (done ? SimTKOpenMMCommon::DefaultReturn : SimTKOpenMMCommon::ErrorReturn);
}
/**---------------------------------------------------------------------------------------
Apply constraint algorithm to velocities.
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param velocities atom velocities
@param inverseMasses 1/mass
@return SimTKOpenMMCommon::DefaultReturn if converge; else
return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int ReferenceShakeAlgorithm::applyToVelocities(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses) {
throw OpenMM::OpenMMException("applyToVelocities is not implemented");
}
/**---------------------------------------------------------------------------------------
Report any violated constriants
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param message report
@return number of violated constraints
--------------------------------------------------------------------------------------- */
int ReferenceShakeAlgorithm::reportShake( int numberOfAtoms, vector<RealVec>& atomCoordinates,
std::stringstream& message ){
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceShakeAlgorithm::reportShake";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM oneM = -1.0;
static const RealOpenMM half = 0.5;
// ---------------------------------------------------------------------------------------
int numberOfConstraints = getNumberOfConstraints();
// loop over constraints calculating distance and comparing to
// expected distance -- report any contraints that are violated
int numberConverged = 0;
RealOpenMM tolerance = getTolerance();
for( int ii = 0; ii < _numberOfConstraints; ii++ ){
int atomI = _atomIndices[ii][0];
int atomJ = _atomIndices[ii][1];
RealOpenMM rp2 = zero;
for( int jj = 0; jj < 3; jj++ ){
rp2 += (atomCoordinates[atomI][jj] - atomCoordinates[atomJ][jj])*(atomCoordinates[atomI][jj] - atomCoordinates[atomJ][jj]);
}
RealOpenMM diff = FABS( rp2 - (_distance[ii]*_distance[ii]) );
if( diff > tolerance ){
message << ii << " constraint violated: " << atomI << " " << atomJ << "] d=" << SQRT( rp2 ) << " " << rp2 << " d0=" << _distance[ii];
message << " diff=" << diff;
message << " [" << atomCoordinates[atomI][0] << " " << atomCoordinates[atomI][1] << " " << atomCoordinates[atomI][2] << "] ";
message << " [" << atomCoordinates[atomJ][0] << " " << atomCoordinates[atomJ][1] << " " << atomCoordinates[atomJ][2] << "] ";
message << "\n";
} else {
numberConverged++;
}
}
return (numberOfConstraints-numberConverged);
}
platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp
View file @ e68df534
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -234,9 +234,8 @@ void ReferenceStochasticDynamics::update(const OpenMM::System& system, vector<Re
updatePart2( numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses, xPrime );
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if( referenceConstraintAlgorithm ){
referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime, inverseMasses );
}
if (referenceConstraintAlgorithm)
referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses);
// copy xPrime -> atomCoordinates
......
platforms/reference/src/SimTKReference/ReferenceVariableStochasticDynamics.cpp
View file @ e68df534
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -278,10 +278,8 @@ void ReferenceVariableStochasticDynamics::update(const OpenMM::System& system, v
updatePart2( numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses, xPrime );
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if( referenceConstraintAlgorithm ){
referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime,
inverseMasses );
}
if (referenceConstraintAlgorithm)
referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses);
// copy xPrime -> atomCoordinates
......
platforms/reference/src/SimTKReference/ReferenceVariableVerletDynamics.cpp
View file @ e68df534
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* Contributors: Peter Eastman, Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -155,7 +155,7 @@ void ReferenceVariableVerletDynamics::update(const OpenMM::System& system, vecto
}
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if (referenceConstraintAlgorithm)
referenceConstraintAlgorithm->apply(numberOfAtoms, atomCoordinates, xPrime, inverseMasses);
referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses);
// Update the positions and velocities.
......
platforms/reference/src/SimTKReference/ReferenceVerletDynamics.cpp
View file @ e68df534
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* Contributors: Peter Eastman, Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -130,8 +130,8 @@ void ReferenceVerletDynamics::update(const OpenMM::System& system, vector<RealVe
}
}
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if( referenceConstraintAlgorithm )
referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime, inverseMasses );
if (referenceConstraintAlgorithm)
referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses);
// Update the positions and velocities.
......
plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp
View file @ e68df534
......@@ -77,7 +77,7 @@ static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>&
if (constraints != NULL) {
constraints->setTolerance(1e-4);
constraints->applyToVelocities(numParticles, posData, shiftedVel, inverseMasses);
constraints->applyToVelocities(posData, shiftedVel, inverseMasses);
}
// Compute the kinetic energy.
......@@ -331,7 +331,7 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
// Apply constraints.
if (constraints != NULL)
constraints->apply(numParticles, pos, xPrime, particleInvMass);
constraints->apply(pos, xPrime, particleInvMass);
// Record the constrained positions and velocities.
......@@ -425,7 +425,7 @@ void ReferenceIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const D
// Apply constraints.
if (constraints != NULL)
constraints->apply(numParticles, pos, xPrime, particleInvMass);
constraints->apply(pos, xPrime, particleInvMass);
// Record the constrained positions and velocities.
......
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