Skip to content

GitLab

Commit bb8c2c6f authored by peastman's avatar peastman
Browse files

Created reference implementation of SETTLE

parent 622cae17
Hide whitespace changes
Inline Side-by-side
Showing with 787 additions and 70 deletions
platforms/reference/include/ReferenceConstraints.h 0 → 100644
View file @ bb8c2c6f
#ifndef OPENMM_REFERENCECONSTRAINTS_H_
#define OPENMM_REFERENCECONSTRAINTS_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ReferenceConstraintAlgorithm.h"
#include "ReferenceCCMAAlgorithm.h"
#include "ReferenceSETTLEAlgorithm.h"
#include "openmm/System.h"
namespace OpenMM {
/**
* This class uses multiple algorithms to apply constraints as efficiently as possible. It identifies clusters
* of three atoms that can be handled by SETTLE, and creates a ReferenceSETTLEAlgorithm object to handle them.
* It then creates a ReferenceCCMAAlgorithm object to handle any remaining constraints.
*/
class OPENMM_EXPORT ReferenceConstraints : public ReferenceConstraintAlgorithm {
public:
ReferenceConstraints(const System& system, RealOpenMM tolerance);
~ReferenceConstraints();
/**
* Get the constraint tolerance.
*/
RealOpenMM getTolerance() const;
/**
* Set the constraint tolerance.
*/
void setTolerance(RealOpenMM tolerance);
/**
* 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);
/**
* 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);
private:
RealOpenMM tolerance;
ReferenceCCMAAlgorithm* ccma;
ReferenceSETTLEAlgorithm* settle;
};
} // namespace OpenMM
#endif /*OPENMM_REFERENCECONSTRAINTS_H_*/
platforms/reference/include/ReferenceSETTLEAlgorithm.h 0 → 100644
View file @ bb8c2c6f
#ifndef OPENMM_REFERENCESETTLEALGORITHM_H_
#define OPENMM_REFERENCESETTLEALGORITHM_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ReferenceConstraintAlgorithm.h"
#include <vector>
namespace OpenMM {
/**
* This implements the SETTLE algorithm for constraining water molecules.
*/
class OPENMM_EXPORT ReferenceSETTLEAlgorithm : public ReferenceConstraintAlgorithm {
public:
ReferenceSETTLEAlgorithm(const std::vector<int>& atom1, const std::vector<int>& atom2, const std::vector<int>& atom3,
const std::vector<RealOpenMM>& distance1, const std::vector<RealOpenMM>& distance2, std::vector<RealOpenMM>& masses, RealOpenMM tolerance);
/**
* Get the constraint tolerance.
*/
RealOpenMM getTolerance() const;
/**
* Set the constraint tolerance.
*/
void setTolerance(RealOpenMM tolerance);
/**
* 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);
/**
* 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);
private:
RealOpenMM tolerance;
std::vector<int> atom1;
std::vector<int> atom2;
std::vector<int> atom3;
std::vector<RealOpenMM> distance1;
std::vector<RealOpenMM> distance2;
std::vector<RealOpenMM> masses;
};
} // namespace OpenMM
#endif /*OPENMM_REFERENCESETTLEALGORITHM_H_*/
platforms/reference/src/ReferenceKernels.cpp
View file @ bb8c2c6f
...@@ -38,6 +38,7 @@ ...@@ -38,6 +38,7 @@
#include "ReferenceBrownianDynamics.h" #include "ReferenceBrownianDynamics.h"
#include "ReferenceCCMAAlgorithm.h" #include "ReferenceCCMAAlgorithm.h"
#include "ReferenceCMAPTorsionIxn.h" #include "ReferenceCMAPTorsionIxn.h"
#include "ReferenceConstraints.h"
#include "ReferenceCustomAngleIxn.h" #include "ReferenceCustomAngleIxn.h"
#include "ReferenceCustomBondIxn.h" #include "ReferenceCustomBondIxn.h"
#include "ReferenceCustomCompoundBondIxn.h" #include "ReferenceCustomCompoundBondIxn.h"
...@@ -132,20 +133,6 @@ static RealVec& extractBoxSize(ContextImpl& context) { ...@@ -132,20 +133,6 @@ static RealVec& extractBoxSize(ContextImpl& context) {
return *(RealVec*) data->periodicBoxSize; return *(RealVec*) data->periodicBoxSize;
} }
static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorithm::AngleInfo>& angles) {
for (int i = 0; i < system.getNumForces(); i++) {
const HarmonicAngleForce* force = dynamic_cast<const HarmonicAngleForce*>(&system.getForce(i));
if (force != NULL) {
for (int j = 0; j < force->getNumAngles(); j++) {
int atom1, atom2, atom3;
double angle, k;
force->getAngleParameters(j, atom1, atom2, atom3, angle, k);
angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, (RealOpenMM)angle));
}
}
}
}
/** /**
* Compute the kinetic energy of the system, possibly shifting the velocities in time to account * Compute the kinetic energy of the system, possibly shifting the velocities in time to account
* for a leapfrog integrator. * for a leapfrog integrator.
...@@ -333,11 +320,8 @@ ReferenceApplyConstraintsKernel::~ReferenceApplyConstraintsKernel() { ...@@ -333,11 +320,8 @@ ReferenceApplyConstraintsKernel::~ReferenceApplyConstraintsKernel() {
} }
void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) { void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
if (constraints == NULL) { if (constraints == NULL)
vector<ReferenceCCMAAlgorithm::AngleInfo> angles; constraints = new ReferenceConstraints(context.getSystem(), (RealOpenMM) tol);
findAnglesForCCMA(context.getSystem(), angles);
constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, tol);
}
vector<RealVec>& positions = extractPositions(context); vector<RealVec>& positions = extractPositions(context);
constraints->setTolerance(tol); constraints->setTolerance(tol);
constraints->apply(data.numParticles, positions, positions, inverseMasses); constraints->apply(data.numParticles, positions, positions, inverseMasses);
...@@ -345,11 +329,8 @@ void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) { ...@@ -345,11 +329,8 @@ void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
} }
void ReferenceApplyConstraintsKernel::applyToVelocities(ContextImpl& context, double tol) { void ReferenceApplyConstraintsKernel::applyToVelocities(ContextImpl& context, double tol) {
if (constraints == NULL) { if (constraints == NULL)
vector<ReferenceCCMAAlgorithm::AngleInfo> angles; constraints = new ReferenceConstraints(context.getSystem(), (RealOpenMM) tol);
findAnglesForCCMA(context.getSystem(), angles);
constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, tol);
}
vector<RealVec>& positions = extractPositions(context); vector<RealVec>& positions = extractPositions(context);
vector<RealVec>& velocities = extractVelocities(context); vector<RealVec>& velocities = extractVelocities(context);
constraints->setTolerance(tol); constraints->setTolerance(tol);
...@@ -1726,9 +1707,7 @@ void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const ...@@ -1726,9 +1707,7 @@ void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const
} }
} }
numConstraints = constraintIndices.size(); numConstraints = constraintIndices.size();
vector<ReferenceCCMAAlgorithm::AngleInfo> angles; constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
findAnglesForCCMA(system, angles);
constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
} }
void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) { void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
...@@ -1780,9 +1759,7 @@ void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, cons ...@@ -1780,9 +1759,7 @@ void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, cons
} }
numConstraints = constraintIndices.size(); numConstraints = constraintIndices.size();
SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed()); SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
vector<ReferenceCCMAAlgorithm::AngleInfo> angles; constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
findAnglesForCCMA(system, angles);
constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
} }
void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) { void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
...@@ -1843,9 +1820,7 @@ void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, cons ...@@ -1843,9 +1820,7 @@ void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, cons
} }
numConstraints = constraintIndices.size(); numConstraints = constraintIndices.size();
SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed()); SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
vector<ReferenceCCMAAlgorithm::AngleInfo> angles; constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
findAnglesForCCMA(system, angles);
constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
} }
void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) { void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
...@@ -1905,9 +1880,7 @@ void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& syst ...@@ -1905,9 +1880,7 @@ void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& syst
} }
numConstraints = constraintIndices.size(); numConstraints = constraintIndices.size();
SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed()); SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
vector<ReferenceCCMAAlgorithm::AngleInfo> angles; constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
findAnglesForCCMA(system, angles);
constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
} }
double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) { double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
...@@ -1967,9 +1940,7 @@ void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system ...@@ -1967,9 +1940,7 @@ void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system
} }
} }
numConstraints = constraintIndices.size(); numConstraints = constraintIndices.size();
vector<ReferenceCCMAAlgorithm::AngleInfo> angles; constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
findAnglesForCCMA(system, angles);
constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
} }
double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) { double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
...@@ -2032,9 +2003,7 @@ void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const ...@@ -2032,9 +2003,7 @@ void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const
dynamics = new ReferenceCustomDynamics(system.getNumParticles(), integrator); dynamics = new ReferenceCustomDynamics(system.getNumParticles(), integrator);
SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed()); SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
vector<ReferenceCCMAAlgorithm::AngleInfo> angles; constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
findAnglesForCCMA(system, angles);
constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
dynamics->setReferenceConstraintAlgorithm(constraints); dynamics->setReferenceConstraintAlgorithm(constraints);
} }
......
platforms/reference/src/SimTKReference/ReferenceConstraints.cpp 0 → 100644
View file @ bb8c2c6f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ReferenceConstraints.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/OpenMMException.h"
#include <map>
#include <utility>
#include <vector>
using namespace OpenMM;
using namespace std;
ReferenceConstraints::ReferenceConstraints(const System& system, RealOpenMM tolerance) : ccma(NULL), settle(NULL) {
this->tolerance = tolerance;
int numParticles = system.getNumParticles();
vector<RealOpenMM> masses(numParticles);
for (int i = 0; i < numParticles; ++i)
masses[i] = (RealOpenMM) system.getParticleMass(i);
// Record the set of constraints and how many constraints each atom is involved in.
vector<int> atom1;
vector<int> atom2;
vector<double> distance;
vector<int> constraintCount(numParticles, 0);
for (int i = 0; i < system.getNumConstraints(); i++) {
int p1, p2;
double d;
system.getConstraintParameters(i, p1, p2, d);
if (masses[p1] != 0 || masses[p2] != 0) {
atom1.push_back(p1);
atom2.push_back(p2);
distance.push_back(d);
constraintCount[p1]++;
constraintCount[p2]++;
}
}
// Identify clusters of three atoms that can be treated with SETTLE. First, for every
// atom that might be part of such a cluster, make a list of the two other atoms it is
// connected to.
vector<map<int, float> > settleConstraints(numParticles);
for (int i = 0; i < (int)atom1.size(); i++) {
if (constraintCount[atom1[i]] == 2 && constraintCount[atom2[i]] == 2) {
settleConstraints[atom1[i]][atom2[i]] = (float) distance[i];
settleConstraints[atom2[i]][atom1[i]] = (float) distance[i];
}
}
// Now remove the ones that don't actually form closed loops of three atoms.
vector<int> settleClusters;
for (int i = 0; i < (int)settleConstraints.size(); i++) {
if (settleConstraints[i].size() == 2) {
int partner1 = settleConstraints[i].begin()->first;
int partner2 = (++settleConstraints[i].begin())->first;
if (settleConstraints[partner1].size() != 2 || settleConstraints[partner2].size() != 2 ||
settleConstraints[partner1].find(partner2) == settleConstraints[partner1].end())
settleConstraints[i].clear();
else if (i < partner1 && i < partner2)
settleClusters.push_back(i);
}
else
settleConstraints[i].clear();
}
// Record the SETTLE clusters.
vector<bool> isSettleAtom(numParticles, false);
int numSETTLE = settleClusters.size();
if (numSETTLE > 0) {
vector<int> atom1(numSETTLE);
vector<int> atom2(numSETTLE);
vector<int> atom3(numSETTLE);
vector<RealOpenMM> distance1(numSETTLE);
vector<RealOpenMM> distance2(numSETTLE);
for (int i = 0; i < numSETTLE; i++) {
int p1 = settleClusters[i];
int p2 = settleConstraints[p1].begin()->first;
int p3 = (++settleConstraints[p1].begin())->first;
float dist12 = settleConstraints[p1].find(p2)->second;
float dist13 = settleConstraints[p1].find(p3)->second;
float dist23 = settleConstraints[p2].find(p3)->second;
if (dist12 == dist13) {
// p1 is the central atom
atom1[i] = p1;
atom2[i] = p2;
atom3[i] = p3;
distance1[i] = dist12;
distance2[i] = dist23;
}
else if (dist12 == dist23) {
// p2 is the central atom
atom1[i] = p2;
atom2[i] = p1;
atom3[i] = p3;
distance1[i] = dist12;
distance2[i] = dist13;
}
else if (dist13 == dist23) {
// p3 is the central atom
atom1[i] = p3;
atom2[i] = p1;
atom3[i] = p2;
distance1[i] = dist13;
distance2[i] = dist12;
}
else
throw OpenMMException("Two of the three distances constrained with SETTLE must be the same.");
isSettleAtom[p1] = true;
isSettleAtom[p2] = true;
isSettleAtom[p3] = true;
}
settle = new ReferenceSETTLEAlgorithm(atom1, atom2, atom3, distance1, distance2, masses, tolerance);
}
// All other constraints are handled with CCMA.
vector<int> ccmaConstraints;
for (unsigned i = 0; i < atom1.size(); i++)
if (!isSettleAtom[atom1[i]])
ccmaConstraints.push_back(i);
int numCCMA = ccmaConstraints.size();
if (numCCMA > 0) {
// Record particles and distances for CCMA.
vector<pair<int, int> > ccmaIndices(numCCMA);
vector<RealOpenMM> ccmaDistance(numCCMA);
for (int i = 0; i < numCCMA; i++) {
int index = ccmaConstraints[i];
ccmaIndices[i] = make_pair(atom1[index], atom2[index]);
ccmaDistance[i] = distance[index];
}
// Look up angles for CCMA.
vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
for (int i = 0; i < system.getNumForces(); i++) {
const HarmonicAngleForce* force = dynamic_cast<const HarmonicAngleForce*>(&system.getForce(i));
if (force != NULL) {
for (int j = 0; j < force->getNumAngles(); j++) {
int atom1, atom2, atom3;
double angle, k;
force->getAngleParameters(j, atom1, atom2, atom3, angle, k);
angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, (RealOpenMM) angle));
}
}
}
// Create the CCMA object.
ccma = new ReferenceCCMAAlgorithm(numParticles, numCCMA, ccmaIndices, ccmaDistance, masses, angles, tolerance);
}
}
ReferenceConstraints::~ReferenceConstraints() {
if (ccma != NULL)
delete ccma;
if (settle != NULL)
delete settle;
}
RealOpenMM ReferenceConstraints::getTolerance() const {
return tolerance;
}
void ReferenceConstraints::setTolerance(RealOpenMM tolerance) {
this->tolerance = tolerance;
if (ccma != NULL)
ccma->setTolerance(tolerance);
if (settle != NULL)
settle->setTolerance(tolerance);
}
int ReferenceConstraints::apply(int numberOfAtoms, vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses) {
if (ccma != NULL)
ccma->apply(numberOfAtoms, atomCoordinates, atomCoordinatesP, inverseMasses);
if (settle != NULL)
settle->apply(numberOfAtoms, atomCoordinates, atomCoordinatesP, inverseMasses);
return SimTKOpenMMCommon::DefaultReturn;
}
int ReferenceConstraints::applyToVelocities(int numberOfAtoms, vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& velocities, vector<RealOpenMM>& inverseMasses) {
if (ccma != NULL)
ccma->applyToVelocities(numberOfAtoms, atomCoordinates, velocities, inverseMasses);
if (settle != NULL)
settle->applyToVelocities(numberOfAtoms, atomCoordinates, velocities, inverseMasses);
return SimTKOpenMMCommon::DefaultReturn;
}
platforms/reference/src/SimTKReference/ReferenceSETTLEAlgorithm.cpp 0 → 100644
View file @ bb8c2c6f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2013 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ReferenceSETTLEAlgorithm.h"
using namespace OpenMM;
using namespace std;
ReferenceSETTLEAlgorithm::ReferenceSETTLEAlgorithm(const vector<int>& atom1, const vector<int>& atom2, const vector<int>& atom3,
const vector<RealOpenMM>& distance1, const vector<RealOpenMM>& distance2, vector<RealOpenMM>& masses, RealOpenMM tolerance) :
atom1(atom1), atom2(atom2), atom3(atom3), distance1(distance1), distance2(distance2), masses(masses), tolerance(tolerance) {
}
RealOpenMM ReferenceSETTLEAlgorithm::getTolerance() const {
return tolerance;
}
void ReferenceSETTLEAlgorithm::setTolerance(RealOpenMM tolerance) {
this->tolerance = tolerance;
}
int ReferenceSETTLEAlgorithm::apply(int numberOfAtoms, 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;
RealVec apos1 = atomCoordinates[atom2[index]];
RealVec xp1 = atomCoordinatesP[atom2[index]]-apos1;
RealVec apos2 = atomCoordinates[atom3[index]];
RealVec xp2 = atomCoordinatesP[atom3[index]]-apos2;
RealOpenMM m0 = masses[atom1[index]];
RealOpenMM m1 = masses[atom2[index]];
RealOpenMM m2 = masses[atom3[index]];
// Apply the SETTLE algorithm.
RealOpenMM xb0 = apos1[0]-apos0[0];
RealOpenMM yb0 = apos1[1]-apos0[1];
RealOpenMM zb0 = apos1[2]-apos0[2];
RealOpenMM xc0 = apos2[0]-apos0[0];
RealOpenMM yc0 = apos2[1]-apos0[1];
RealOpenMM zc0 = apos2[2]-apos0[2];
RealOpenMM invTotalMass = 1/(m0+m1+m2);
RealOpenMM xcom = (xp0[0]*m0 + (xb0+xp1[0])*m1 + (xc0+xp2[0])*m2) * invTotalMass;
RealOpenMM ycom = (xp0[1]*m0 + (yb0+xp1[1])*m1 + (yc0+xp2[1])*m2) * invTotalMass;
RealOpenMM zcom = (xp0[2]*m0 + (zb0+xp1[2])*m1 + (zc0+xp2[2])*m2) * invTotalMass;
RealOpenMM xa1 = xp0[0] - xcom;
RealOpenMM ya1 = xp0[1] - ycom;
RealOpenMM za1 = xp0[2] - zcom;
RealOpenMM xb1 = xb0 + xp1[0] - xcom;
RealOpenMM yb1 = yb0 + xp1[1] - ycom;
RealOpenMM zb1 = zb0 + xp1[2] - zcom;
RealOpenMM xc1 = xc0 + xp2[0] - xcom;
RealOpenMM yc1 = yc0 + xp2[1] - ycom;
RealOpenMM zc1 = zc0 + xp2[2] - zcom;
RealOpenMM xaksZd = yb0*zc0 - zb0*yc0;
RealOpenMM yaksZd = zb0*xc0 - xb0*zc0;
RealOpenMM zaksZd = xb0*yc0 - yb0*xc0;
RealOpenMM xaksXd = ya1*zaksZd - za1*yaksZd;
RealOpenMM yaksXd = za1*xaksZd - xa1*zaksZd;
RealOpenMM zaksXd = xa1*yaksZd - ya1*xaksZd;
RealOpenMM xaksYd = yaksZd*zaksXd - zaksZd*yaksXd;
RealOpenMM yaksYd = zaksZd*xaksXd - xaksZd*zaksXd;
RealOpenMM zaksYd = xaksZd*yaksXd - yaksZd*xaksXd;
RealOpenMM axlng = sqrt(xaksXd*xaksXd + yaksXd*yaksXd + zaksXd*zaksXd);
RealOpenMM aylng = sqrt(xaksYd*xaksYd + yaksYd*yaksYd + zaksYd*zaksYd);
RealOpenMM azlng = sqrt(xaksZd*xaksZd + yaksZd*yaksZd + zaksZd*zaksZd);
RealOpenMM trns11 = xaksXd / axlng;
RealOpenMM trns21 = yaksXd / axlng;
RealOpenMM trns31 = zaksXd / axlng;
RealOpenMM trns12 = xaksYd / aylng;
RealOpenMM trns22 = yaksYd / aylng;
RealOpenMM trns32 = zaksYd / aylng;
RealOpenMM trns13 = xaksZd / azlng;
RealOpenMM trns23 = yaksZd / azlng;
RealOpenMM trns33 = zaksZd / azlng;
RealOpenMM xb0d = trns11*xb0 + trns21*yb0 + trns31*zb0;
RealOpenMM yb0d = trns12*xb0 + trns22*yb0 + trns32*zb0;
RealOpenMM xc0d = trns11*xc0 + trns21*yc0 + trns31*zc0;
RealOpenMM yc0d = trns12*xc0 + trns22*yc0 + trns32*zc0;
RealOpenMM za1d = trns13*xa1 + trns23*ya1 + trns33*za1;
RealOpenMM xb1d = trns11*xb1 + trns21*yb1 + trns31*zb1;
RealOpenMM yb1d = trns12*xb1 + trns22*yb1 + trns32*zb1;
RealOpenMM zb1d = trns13*xb1 + trns23*yb1 + trns33*zb1;
RealOpenMM xc1d = trns11*xc1 + trns21*yc1 + trns31*zc1;
RealOpenMM yc1d = trns12*xc1 + trns22*yc1 + trns32*zc1;
RealOpenMM zc1d = trns13*xc1 + trns23*yc1 + trns33*zc1;
// --- Step2 A2' ---
RealOpenMM rc = 0.5*distance2[index];
RealOpenMM rb = sqrt(distance1[index]*distance1[index]-rc*rc);
RealOpenMM ra = rb*(m1+m2)*invTotalMass;
rb -= ra;
RealOpenMM sinphi = za1d / ra;
RealOpenMM cosphi = sqrt(1 - sinphi*sinphi);
RealOpenMM sinpsi = (zb1d - zc1d) / (2*rc*cosphi);
RealOpenMM cospsi = sqrt(1 - sinpsi*sinpsi);
RealOpenMM ya2d = ra*cosphi;
RealOpenMM xb2d = - rc*cospsi;
RealOpenMM yb2d = - rb*cosphi - rc*sinpsi*sinphi;
RealOpenMM yc2d = - rb*cosphi + rc*sinpsi*sinphi;
RealOpenMM xb2d2 = xb2d*xb2d;
RealOpenMM hh2 = 4.0f*xb2d2 + (yb2d-yc2d)*(yb2d-yc2d) + (zb1d-zc1d)*(zb1d-zc1d);
RealOpenMM deltx = 2.0f*xb2d + sqrt(4.0f*xb2d2 - hh2 + distance2[index]*distance2[index]);
xb2d -= deltx*0.5;
// --- Step3 al,be,ga ---
RealOpenMM alpha = (xb2d*(xb0d-xc0d) + yb0d*yb2d + yc0d*yc2d);
RealOpenMM beta = (xb2d*(yc0d-yb0d) + xb0d*yb2d + xc0d*yc2d);
RealOpenMM gamma = xb0d*yb1d - xb1d*yb0d + xc0d*yc1d - xc1d*yc0d;
RealOpenMM al2be2 = alpha*alpha + beta*beta;
RealOpenMM sintheta = (alpha*gamma - beta*sqrt(al2be2 - gamma*gamma)) / al2be2;
// --- Step4 A3' ---
RealOpenMM costheta = sqrt(1 - sintheta*sintheta);
RealOpenMM xa3d = - ya2d*sintheta;
RealOpenMM ya3d = ya2d*costheta;
RealOpenMM za3d = za1d;
RealOpenMM xb3d = xb2d*costheta - yb2d*sintheta;
RealOpenMM yb3d = xb2d*sintheta + yb2d*costheta;
RealOpenMM zb3d = zb1d;
RealOpenMM xc3d = - xb2d*costheta - yc2d*sintheta;
RealOpenMM yc3d = - xb2d*sintheta + yc2d*costheta;
RealOpenMM zc3d = zc1d;
// --- Step5 A3 ---
RealOpenMM xa3 = trns11*xa3d + trns12*ya3d + trns13*za3d;
RealOpenMM ya3 = trns21*xa3d + trns22*ya3d + trns23*za3d;
RealOpenMM za3 = trns31*xa3d + trns32*ya3d + trns33*za3d;
RealOpenMM xb3 = trns11*xb3d + trns12*yb3d + trns13*zb3d;
RealOpenMM yb3 = trns21*xb3d + trns22*yb3d + trns23*zb3d;
RealOpenMM zb3 = trns31*xb3d + trns32*yb3d + trns33*zb3d;
RealOpenMM xc3 = trns11*xc3d + trns12*yc3d + trns13*zc3d;
RealOpenMM yc3 = trns21*xc3d + trns22*yc3d + trns23*zc3d;
RealOpenMM zc3 = trns31*xc3d + trns32*yc3d + trns33*zc3d;
xp0[0] = xcom + xa3;
xp0[1] = ycom + ya3;
xp0[2] = zcom + za3;
xp1[0] = xcom + xb3 - xb0;
xp1[1] = ycom + yb3 - yb0;
xp1[2] = zcom + zb3 - zb0;
xp2[0] = xcom + xc3 - xc0;
xp2[1] = ycom + yc3 - yc0;
xp2[2] = zcom + zc3 - zc0;
// Record the new positions.
atomCoordinatesP[atom1[index]] = xp0+apos0;
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) {
for (int index = 0; index < (int) atom1.size(); ++index) {
RealVec apos0 = atomCoordinates[atom1[index]];
RealVec apos1 = atomCoordinates[atom2[index]];
RealVec apos2 = atomCoordinates[atom3[index]];
RealVec v0 = velocities[atom1[index]];
RealVec v1 = velocities[atom2[index]];
RealVec v2 = velocities[atom3[index]];
// Compute intermediate quantities: the atom masses, the bond directions, the relative velocities,
// and the angle cosines and sines.
RealOpenMM mA = masses[atom1[index]];
RealOpenMM mB = masses[atom2[index]];
RealOpenMM mC = masses[atom3[index]];
RealVec eAB = apos1-apos0;
RealVec eBC = apos2-apos1;
RealVec eCA = apos0-apos2;
eAB /= sqrt(eAB[0]*eAB[0] + eAB[1]*eAB[1] + eAB[2]*eAB[2]);
eBC /= sqrt(eBC[0]*eBC[0] + eBC[1]*eBC[1] + eBC[2]*eBC[2]);
eCA /= sqrt(eCA[0]*eCA[0] + eCA[1]*eCA[1] + eCA[2]*eCA[2]);
RealOpenMM vAB = (v1[0]-v0[0])*eAB[0] + (v1[1]-v0[1])*eAB[1] + (v1[2]-v0[2])*eAB[2];
RealOpenMM vBC = (v2[0]-v1[0])*eBC[0] + (v2[1]-v1[1])*eBC[1] + (v2[2]-v1[2])*eBC[2];
RealOpenMM vCA = (v0[0]-v2[0])*eCA[0] + (v0[1]-v2[1])*eCA[1] + (v0[2]-v2[2])*eCA[2];
RealOpenMM cA = -(eAB[0]*eCA[0] + eAB[1]*eCA[1] + eAB[2]*eCA[2]);
RealOpenMM cB = -(eAB[0]*eBC[0] + eAB[1]*eBC[1] + eAB[2]*eBC[2]);
RealOpenMM cC = -(eBC[0]*eCA[0] + eBC[1]*eCA[1] + eBC[2]*eCA[2]);
RealOpenMM s2A = 1-cA*cA;
RealOpenMM s2B = 1-cB*cB;
RealOpenMM s2C = 1-cC*cC;
// Solve the equations. These are different from those in the SETTLE paper (JCC 13(8), pp. 952-962, 1992), because
// in going from equations B1 to B2, they make the assumption that mB=mC (but don't bother to mention they're
// making that assumption). We allow all three atoms to have different masses.
RealOpenMM mABCinv = 1/(mA*mB*mC);
RealOpenMM denom = (((s2A*mB+s2B*mA)*mC+(s2A*mB*mB+2*(cA*cB*cC+1)*mA*mB+s2B*mA*mA))*mC+s2C*mA*mB*(mA+mB))*mABCinv;
RealOpenMM tab = ((cB*cC*mA-cA*mB-cA*mC)*vCA + (cA*cC*mB-cB*mC-cB*mA)*vBC + (s2C*mA*mA*mB*mB*mABCinv+(mA+mB+mC))*vAB)/denom;
RealOpenMM tbc = ((cA*cB*mC-cC*mB-cC*mA)*vCA + (s2A*mB*mB*mC*mC*mABCinv+(mA+mB+mC))*vBC + (cA*cC*mB-cB*mA-cB*mC)*vAB)/denom;
RealOpenMM tca = ((s2B*mA*mA*mC*mC*mABCinv+(mA+mB+mC))*vCA + (cA*cB*mC-cC*mB-cC*mA)*vBC + (cB*cC*mA-cA*mB-cA*mC)*vAB)/denom;
v0 += (eAB*tab - eCA*tca)*inverseMasses[atom1[index]];
v1 += (eBC*tbc - eAB*tab)*inverseMasses[atom2[index]];
v2 += (eCA*tca - eBC*tbc)*inverseMasses[atom3[index]];
velocities[atom1[index]] = v0;
velocities[atom2[index]] = v1;
velocities[atom3[index]] = v2;
}
return SimTKOpenMMCommon::DefaultReturn;
}
platforms/reference/tests/TestReferenceCustomIntegrator.cpp
View file @ bb8c2c6f
...@@ -159,7 +159,7 @@ void testConstraints() { ...@@ -159,7 +159,7 @@ void testConstraints() {
* Test an integrator that applies constraints directly to velocities. * Test an integrator that applies constraints directly to velocities.
*/ */
void testVelocityConstraints() { void testVelocityConstraints() {
const int numParticles = 8; const int numParticles = 10;
ReferencePlatform platform; ReferencePlatform platform;
System system; System system;
CustomIntegrator integrator(0.002); CustomIntegrator integrator(0.002);
...@@ -176,7 +176,21 @@ void testVelocityConstraints() { ...@@ -176,7 +176,21 @@ void testVelocityConstraints() {
system.addParticle(i%2 == 0 ? 5.0 : 10.0); system.addParticle(i%2 == 0 ? 5.0 : 10.0);
forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0); forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
} }
for (int i = 0; i < numParticles-1; ++i)
// Constrain the first three particles with SHAKE.
system.addConstraint(0, 1, 1.0);
system.addConstraint(1, 2, 1.0);
// Constrain the next three with SETTLE.
system.addConstraint(3, 4, 1.0);
system.addConstraint(5, 4, 1.0);
system.addConstraint(3, 5, sqrt(2.0));
// Constraint the rest with CCMA.
for (int i = 6; i < numParticles-1; ++i)
system.addConstraint(i, i+1, 1.0); system.addConstraint(i, i+1, 1.0);
system.addForce(forceField); system.addForce(forceField);
Context context(system, integrator, platform); Context context(system, integrator, platform);
...@@ -196,6 +210,7 @@ void testVelocityConstraints() { ...@@ -196,6 +210,7 @@ void testVelocityConstraints() {
double initialEnergy = 0.0; double initialEnergy = 0.0;
for (int i = 0; i < 1000; ++i) { for (int i = 0; i < 1000; ++i) {
integrator.step(2);
State state = context.getState(State::Positions | State::Velocities | State::Energy); State state = context.getState(State::Positions | State::Velocities | State::Energy);
for (int j = 0; j < system.getNumConstraints(); ++j) { for (int j = 0; j < system.getNumConstraints(); ++j) {
int particle1, particle2; int particle1, particle2;
...@@ -213,11 +228,10 @@ void testVelocityConstraints() { ...@@ -213,11 +228,10 @@ void testVelocityConstraints() {
} }
} }
double energy = state.getKineticEnergy()+state.getPotentialEnergy(); double energy = state.getKineticEnergy()+state.getPotentialEnergy();
if (i == 1) if (i == 0)
initialEnergy = energy; initialEnergy = energy;
else if (i > 1) else if (i > 0)
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01); ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
integrator.step(2);
} }
} }
......
platforms/reference/tests/TestReferenceSettle.cpp 0 → 100644
View file @ bb8c2c6f
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of the SETTLE algorithm.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
void testConstraints() {
const int numMolecules = 10;
const int numParticles = numMolecules*3;
const int numConstraints = numMolecules*3;
const double temp = 100.0;
ReferencePlatform platform;
System system;
LangevinIntegrator integrator(temp, 2.0, 0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numMolecules; ++i) {
system.addParticle(16.0);
system.addParticle(1.0);
system.addParticle(1.0);
forceField->addParticle(-0.82, 0.317, 0.65);
forceField->addParticle(0.41, 1.0, 0.0);
forceField->addParticle(0.41, 1.0, 0.0);
system.addConstraint(i*3, i*3+1, 0.1);
system.addConstraint(i*3, i*3+2, 0.1);
system.addConstraint(i*3+1, i*3+2, 0.163);
}
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numMolecules; ++i) {
positions[i*3] = Vec3((i%4)*0.4, (i/4)*0.4, 0);
positions[i*3+1] = positions[i*3]+Vec3(0.1, 0, 0);
positions[i*3+2] = positions[i*3]+Vec3(-0.03333, 0.09428, 0);
velocities[i*3] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
velocities[i*3+1] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
velocities[i*3+2] = Vec3(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
}
context.setPositions(positions);
context.setVelocities(velocities);
// Simulate it and see whether the constraints remain satisfied.
for (int i = 0; i < 1000; ++i) {
integrator.step(1);
State state = context.getState(State::Positions | State::Forces);
for (int j = 0; j < numConstraints; ++j) {
int particle1, particle2;
double distance;
system.getConstraintParameters(j, particle1, particle2, distance);
Vec3 p1 = state.getPositions()[particle1];
Vec3 p2 = state.getPositions()[particle2];
double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
ASSERT_EQUAL_TOL(distance, dist, 1e-5);
}
}
}
int main(int argc, char* argv[]) {
try {
testConstraints();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp
View file @ bb8c2c6f
...@@ -34,7 +34,7 @@ ...@@ -34,7 +34,7 @@
#include "openmm/OpenMMException.h" #include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h" #include "openmm/internal/ContextImpl.h"
#include "SimTKOpenMMUtilities.h" #include "SimTKOpenMMUtilities.h"
#include "ReferenceCCMAAlgorithm.h" #include "ReferenceConstraints.h"
#include "ReferenceVirtualSites.h" #include "ReferenceVirtualSites.h"
#include <set> #include <set>
...@@ -56,20 +56,6 @@ static vector<RealVec>& extractForces(ContextImpl& context) { ...@@ -56,20 +56,6 @@ static vector<RealVec>& extractForces(ContextImpl& context) {
return *((vector<RealVec>*) data->forces); return *((vector<RealVec>*) data->forces);
} }
static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorithm::AngleInfo>& angles) {
for (int i = 0; i < system.getNumForces(); i++) {
const HarmonicAngleForce* force = dynamic_cast<const HarmonicAngleForce*>(&system.getForce(i));
if (force != NULL) {
for (int j = 0; j < force->getNumAngles(); j++) {
int atom1, atom2, atom3;
double angle, k;
force->getAngleParameters(j, atom1, atom2, atom3, angle, k);
angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, (RealOpenMM)angle));
}
}
}
}
static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& inverseMasses, double timeShift, ReferenceConstraintAlgorithm* constraints) { static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& inverseMasses, double timeShift, ReferenceConstraintAlgorithm* constraints) {
const System& system = context.getSystem(); const System& system = context.getSystem();
int numParticles = system.getNumParticles(); int numParticles = system.getNumParticles();
...@@ -283,10 +269,7 @@ void ReferenceIntegrateDrudeLangevinStepKernel::initialize(const System& system, ...@@ -283,10 +269,7 @@ void ReferenceIntegrateDrudeLangevinStepKernel::initialize(const System& system,
constraintDistances.push_back(distance); constraintDistances.push_back(distance);
} }
} }
int numConstraints = constraintIndices.size(); constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
findAnglesForCCMA(system, angles);
constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, particleMass, angles, (RealOpenMM)integrator.getConstraintTolerance());
} }
} }
...@@ -408,10 +391,7 @@ void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, cons ...@@ -408,10 +391,7 @@ void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, cons
constraintDistances.push_back(distance); constraintDistances.push_back(distance);
} }
} }
int numConstraints = constraintIndices.size(); constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
findAnglesForCCMA(system, angles);
constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, particleMass, angles, (RealOpenMM)integrator.getConstraintTolerance());
} }
// Initialize the energy minimizer. // Initialize the energy minimizer.
......
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