Skip to content

GitLab

Commit 64b39ffd authored by peastman's avatar peastman
Browse files

Merge pull request #164 from peastman/master 

Lots of code cleanup related to constraints in the reference platform
parents bf8db55f e68df534
Show whitespace changes
Inline Side-by-side
Showing with 366 additions and 1628 deletions
platforms/reference/include/ReferenceCCMAAlgorithm.h
View file @ 64b39ffd
/* 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 @ bf8db55f
/* 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 @ 64b39ffd
/* 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 @ 64b39ffd
......@@ -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 @ 64b39ffd
......@@ -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 @ bf8db55f
/* 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 @ 64b39ffd
......@@ -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 @ 64b39ffd
/* 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 @ 64b39ffd
/* 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 @ bf8db55f
/* 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 @ 64b39ffd
......@@ -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 @ 64b39ffd
/* 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 @ 64b39ffd
......@@ -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 @ bf8db55f
/* 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 @ 64b39ffd
/* 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 @ 64b39ffd
/* 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 @ 64b39ffd
/* 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 @ 64b39ffd
/* 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 @ 64b39ffd
......@@ -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