Skip to content

GitLab

Commit fd01aa58 authored by Peter Eastman's avatar Peter Eastman
Browse files

Implemented virtual sites for reference platform

parent cbd45dca
Hide whitespace changes
Inline Side-by-side
Showing with 683 additions and 272 deletions
platforms/reference/src/ReferenceKernels.cpp
View file @ fd01aa58
......@@ -6,7 +6,7 @@
* 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. *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -56,6 +56,7 @@
#include "SimTKReference/ReferenceVariableStochasticDynamics.h"
#include "SimTKReference/ReferenceVariableVerletDynamics.h"
#include "SimTKReference/ReferenceVerletDynamics.h"
#include "SimTKReference/ReferenceVirtualSites.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/Context.h"
#include "openmm/System.h"
......@@ -161,6 +162,8 @@ void ReferenceCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context,
double ReferenceCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy) {
if (!includeForces)
extractForces(context) = savedForces; // Restore the forces so computing the energy doesn't overwrite the forces with incorrect values.
else
ReferenceVirtualSites::distributeForces(context.getSystem(), extractPositions(context), extractForces(context));
return 0.0;
}
......@@ -272,6 +275,7 @@ void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
vector<RealVec>& positions = extractPositions(context);
constraints->setTolerance(tol);
constraints->apply(data.numParticles, positions, positions, inverseMasses);
ReferenceVirtualSites::computePositions(context.getSystem(), positions);
}
ReferenceCalcHarmonicBondForceKernel::~ReferenceCalcHarmonicBondForceKernel() {
......@@ -1261,7 +1265,7 @@ void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const Ver
prevStepSize = stepSize;
}
constraints->setTolerance(integrator.getConstraintTolerance());
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
dynamics->update(context.getSystem(), posData, velData, forceData, masses);
data.time += stepSize;
data.stepCount++;
}
......@@ -1325,7 +1329,7 @@ void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const L
prevStepSize = stepSize;
}
constraints->setTolerance(integrator.getConstraintTolerance());
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
dynamics->update(context.getSystem(), posData, velData, forceData, masses);
data.time += stepSize;
data.stepCount++;
}
......@@ -1388,7 +1392,7 @@ void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const B
prevStepSize = stepSize;
}
constraints->setTolerance(integrator.getConstraintTolerance());
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
dynamics->update(context.getSystem(), posData, velData, forceData, masses);
data.time += stepSize;
data.stepCount++;
}
......@@ -1449,7 +1453,7 @@ void ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context,
}
constraints->setTolerance(integrator.getConstraintTolerance());
RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize);
data.time += dynamics->getDeltaT();
if (dynamics->getDeltaT() == maxStepSize)
data.time = maxTime; // Avoid round-off error
......@@ -1506,7 +1510,7 @@ void ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, c
}
constraints->setTolerance(integrator.getConstraintTolerance());
RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize);
data.time += dynamics->getDeltaT();
if (dynamics->getDeltaT() == maxStepSize)
data.time = maxTime; // Avoid round-off error
......@@ -1676,10 +1680,10 @@ void ReferenceRemoveCMMotionKernel::execute(ContextImpl& context) {
RealOpenMM momentum[] = {0.0, 0.0, 0.0};
RealOpenMM mass = 0.0;
for (size_t i = 0; i < masses.size(); ++i) {
momentum[0] += static_cast<RealOpenMM>( masses[i]*velData[i][0] );
momentum[1] += static_cast<RealOpenMM>( masses[i]*velData[i][1] );
momentum[2] += static_cast<RealOpenMM>( masses[i]*velData[i][2] );
mass += static_cast<RealOpenMM>( masses[i] );
momentum[0] += static_cast<RealOpenMM>(masses[i]*velData[i][0]);
momentum[1] += static_cast<RealOpenMM>(masses[i]*velData[i][1]);
momentum[2] += static_cast<RealOpenMM>(masses[i]*velData[i][2]);
mass += static_cast<RealOpenMM>(masses[i]);
}
// Adjust the particle velocities.
......@@ -1688,8 +1692,10 @@ void ReferenceRemoveCMMotionKernel::execute(ContextImpl& context) {
momentum[1] /= mass;
momentum[2] /= mass;
for (size_t i = 0; i < masses.size(); ++i) {
velData[i][0] -= momentum[0];
velData[i][1] -= momentum[1];
velData[i][2] -= momentum[2];
if (masses[i] != 0.0) {
velData[i][0] -= momentum[0];
velData[i][1] -= momentum[1];
velData[i][2] -= momentum[2];
}
}
}
platforms/reference/src/SimTKReference/ReferenceBrownianDynamics.cpp
View file @ fd01aa58
/* Portions copyright (c) 2006-2008 Stanford University and Simbios.
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Peter Eastman, Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -29,6 +29,7 @@
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceBrownianDynamics.h"
#include "ReferenceVirtualSites.h"
#include <cstdio>
......@@ -114,7 +115,7 @@ RealOpenMM ReferenceBrownianDynamics::getFriction( void ) const {
Update -- driver routine for performing Brownian dynamics update of coordinates
and velocities
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -122,7 +123,7 @@ RealOpenMM ReferenceBrownianDynamics::getFriction( void ) const {
--------------------------------------------------------------------------------------- */
void ReferenceBrownianDynamics::update( int numberOfAtoms, vector<RealVec>& atomCoordinates,
void ReferenceBrownianDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& masses ){
......@@ -137,27 +138,15 @@ void ReferenceBrownianDynamics::update( int numberOfAtoms, vector<RealVec>& atom
// first-time-through initialization
int numberOfAtoms = system.getNumParticles();
if( getTimeStep() == 0 ){
std::stringstream message;
message << methodName;
int errors = 0;
// invert masses
for( int ii = 0; ii < numberOfAtoms; ii++ ){
if( masses[ii] <= zero ){
message << "mass at atom index=" << ii << " (" << masses[ii] << ") is <= 0" << std::endl;
errors++;
} else {
inverseMasses[ii] = one/masses[ii];
}
}
// exit if errors
if( errors ){
SimTKOpenMMLog::printError( message );
if (masses[ii] == zero)
inverseMasses[ii] == zero;
else
inverseMasses[ii] = one/masses[ii];
}
}
......@@ -166,9 +155,10 @@ void ReferenceBrownianDynamics::update( int numberOfAtoms, vector<RealVec>& atom
const RealOpenMM noiseAmplitude = static_cast<RealOpenMM>( sqrt(2.0*BOLTZ*getTemperature()*getDeltaT()/getFriction()) );
const RealOpenMM forceScale = getDeltaT()/getFriction();
for (int i = 0; i < numberOfAtoms; ++i) {
for (int j = 0; j < 3; ++j) {
xPrime[i][j] = atomCoordinates[i][j] + forceScale*inverseMasses[i]*forces[i][j] + noiseAmplitude*SQRT(inverseMasses[i])*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
if (masses[i] != zero)
for (int j = 0; j < 3; ++j) {
xPrime[i][j] = atomCoordinates[i][j] + forceScale*inverseMasses[i]*forces[i][j] + noiseAmplitude*SQRT(inverseMasses[i])*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
}
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if( referenceConstraintAlgorithm )
......@@ -178,12 +168,12 @@ void ReferenceBrownianDynamics::update( int numberOfAtoms, vector<RealVec>& atom
RealOpenMM velocityScale = static_cast<RealOpenMM>( 1.0/getDeltaT() );
for (int i = 0; i < numberOfAtoms; ++i) {
for (int j = 0; j < 3; ++j) {
velocities[i][j] = velocityScale*(xPrime[i][j] - atomCoordinates[i][j]);
atomCoordinates[i][j] = xPrime[i][j];
}
if (masses[i] != zero)
for (int j = 0; j < 3; ++j) {
velocities[i][j] = velocityScale*(xPrime[i][j] - atomCoordinates[i][j]);
atomCoordinates[i][j] = xPrime[i][j];
}
}
ReferenceVirtualSites::computePositions(system, atomCoordinates);
incrementTimeStep();
}
platforms/reference/src/SimTKReference/ReferenceBrownianDynamics.h
View file @ fd01aa58
/* Portions copyright (c) 2006-2008 Stanford University and Simbios.
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -74,7 +74,7 @@ class ReferenceBrownianDynamics : public ReferenceDynamics {
Update
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -82,7 +82,7 @@ class ReferenceBrownianDynamics : public ReferenceDynamics {
--------------------------------------------------------------------------------------- */
void update( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses );
};
......
platforms/reference/src/SimTKReference/ReferenceCustomDynamics.cpp
View file @ fd01aa58
......@@ -25,13 +25,14 @@
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "openmm/OpenMMException.h"
#include "lepton/Parser.h"
#include "ReferenceVirtualSites.h"
#include "ReferenceCustomDynamics.h"
#include "lepton/ParsedExpression.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/ForceImpl.h"
#include "lepton/Operation.h"
#include "lepton/ParsedExpression.h"
#include "lepton/Parser.h"
#include <set>
using namespace std;
......@@ -126,8 +127,12 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
// Build the list of inverse masses.
inverseMasses.resize(numberOfAtoms);
for (int i = 0; i < numberOfAtoms; i++)
inverseMasses[i] = 1.0/masses[i];
for (int i = 0; i < numberOfAtoms; i++) {
if (masses[i] == 0.0)
inverseMasses[i] = 0.0;
else
inverseMasses[i] = 1.0/masses[i];
}
}
// Loop over steps and execute them.
......@@ -188,7 +193,8 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
computePerDof(numberOfAtoms, sumBuffer, atomCoordinates, velocities, forces, masses, globals, perDof, stepExpression[i]);
RealOpenMM sum = 0.0;
for (int j = 0; j < numberOfAtoms; j++)
sum += sumBuffer[j][0]+sumBuffer[j][1]+sumBuffer[j][2];
if (masses[j] != 0.0)
sum += sumBuffer[j][0]+sumBuffer[j][1]+sumBuffer[j][2];
globals[stepVariable[i]] = sum;
break;
}
......@@ -209,6 +215,7 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
if (invalidatesForces[i])
forcesAreValid = false;
}
ReferenceVirtualSites::computePositions(context.getSystem(), atomCoordinates);
incrementTimeStep();
recordChangedParameters(context, globals);
}
......@@ -220,18 +227,20 @@ RealOpenMM ReferenceCustomDynamics::computePerDof(int numberOfAtoms, vector<Real
map<string, RealOpenMM> variables = globals;
for (int i = 0; i < numberOfAtoms; i++) {
variables["m"] = masses[i];
for (int j = 0; j < 3; j++) {
// Compute the expression.
variables["x"] = atomCoordinates[i][j];
variables["v"] = velocities[i][j];
variables["f"] = forces[i][j];
variables["uniform"] = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber();
variables["gaussian"] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
for (int k = 0; k < (int) perDof.size(); k++)
variables[integrator.getPerDofVariableName(k)] = perDof[k][i][j];
results[i][j] = expression.evaluate(variables);
if (masses[i] != 0.0) {
variables["m"] = masses[i];
for (int j = 0; j < 3; j++) {
// Compute the expression.
variables["x"] = atomCoordinates[i][j];
variables["v"] = velocities[i][j];
variables["f"] = forces[i][j];
variables["uniform"] = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber();
variables["gaussian"] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
for (int k = 0; k < (int) perDof.size(); k++)
variables[integrator.getPerDofVariableName(k)] = perDof[k][i][j];
results[i][j] = expression.evaluate(variables);
}
}
}
}
......
platforms/reference/src/SimTKReference/ReferenceDynamics.cpp
View file @ fd01aa58
......@@ -241,10 +241,10 @@ void ReferenceDynamics::setReferenceConstraintAlgorithm( ReferenceConstraintAlgo
/**---------------------------------------------------------------------------------------
Update -- driver routine for performing stochastic dynamics update of coordinates
Update -- driver routine for performing dynamics update of coordinates
and velocities
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -252,7 +252,7 @@ void ReferenceDynamics::setReferenceConstraintAlgorithm( ReferenceConstraintAlgo
--------------------------------------------------------------------------------------- */
void ReferenceDynamics::update( int numberOfAtoms, vector<RealVec>& atomCoordinates,
void ReferenceDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses ){
// ---------------------------------------------------------------------------------------
......@@ -264,49 +264,3 @@ void ReferenceDynamics::update( int numberOfAtoms, vector<RealVec>& atomCoordina
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
Remove total linear momentum
@param numberOfAtoms number of atoms
@param masses masses
@param velocities velocities
--------------------------------------------------------------------------------------- */
void ReferenceDynamics::removeTotalLinearMomentum( int numberOfAtoms, RealOpenMM* masses,
vector<RealVec>& velocities ) const {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceDynamics::removeTotalLinearMomentum";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
RealOpenMM totalMass = zero;
RealOpenMM linearMomentum[3] = { zero, zero, zero };
for( int ii = 0; ii < numberOfAtoms; ii++ ){
totalMass += masses[ii];
linearMomentum[0] += masses[ii]*velocities[ii][0];
linearMomentum[1] += masses[ii]*velocities[ii][1];
linearMomentum[2] += masses[ii]*velocities[ii][2];
}
if( totalMass > zero ){
RealOpenMM totalMassI = one/totalMass;
linearMomentum[0] *= totalMassI;
linearMomentum[1] *= totalMassI;
linearMomentum[2] *= totalMassI;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
velocities[ii][0] -= linearMomentum[0];
velocities[ii][1] -= linearMomentum[1];
velocities[ii][2] -= linearMomentum[2];
}
}
}
platforms/reference/src/SimTKReference/ReferenceDynamics.h
View file @ fd01aa58
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -27,6 +27,7 @@
#include "ReferenceConstraintAlgorithm.h"
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "openmm/System.h"
#include <cstddef>
#include <vector>
......@@ -134,24 +135,12 @@ class OPENMM_EXPORT ReferenceDynamics {
--------------------------------------------------------------------------------------- */
RealOpenMM getTemperature( void ) const;
/**---------------------------------------------------------------------------------------
Remove total linear momentum
@param numberOfAtoms number of atoms
@param masses masses
@param velocities velocities
--------------------------------------------------------------------------------------- */
void removeTotalLinearMomentum( int numberOfAtoms, RealOpenMM* masses, std::vector<OpenMM::RealVec>& velocities ) const;
/**---------------------------------------------------------------------------------------
Update
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -159,8 +148,8 @@ class OPENMM_EXPORT ReferenceDynamics {
--------------------------------------------------------------------------------------- */
virtual void update( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses );
virtual void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses);
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp
View file @ fd01aa58
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -29,6 +29,7 @@
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceStochasticDynamics.h"
#include "ReferenceVirtualSites.h"
#include <cstdio>
......@@ -144,10 +145,12 @@ void ReferenceStochasticDynamics::updatePart1( int numberOfAtoms, vector<RealVec
const RealOpenMM noisescale = SQRT(2*kT/tau)*SQRT(0.5*(1-vscale*vscale)*tau);
for (int ii = 0; ii < numberOfAtoms; ii++) {
RealOpenMM sqrtInvMass = SQRT(inverseMasses[ii]);
for (int jj = 0; jj < 3; jj++) {
velocities[ii][jj] = vscale*velocities[ii][jj] + fscale*inverseMasses[ii]*forces[ii][jj] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
if (inverseMasses[ii] != 0.0) {
RealOpenMM sqrtInvMass = SQRT(inverseMasses[ii]);
for (int jj = 0; jj < 3; jj++) {
velocities[ii][jj] = vscale*velocities[ii][jj] + fscale*inverseMasses[ii]*forces[ii][jj] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
}
}
}
......@@ -177,9 +180,9 @@ void ReferenceStochasticDynamics::updatePart2( int numberOfAtoms, vector<RealVec
// perform second update
for (int ii = 0; ii < numberOfAtoms; ii++) {
for (int jj = 0; jj < 3; jj++) {
xPrime[ii][jj] = atomCoordinates[ii][jj]+getDeltaT()*velocities[ii][jj];
}
if (inverseMasses[ii] != 0.0)
for (int jj = 0; jj < 3; jj++)
xPrime[ii][jj] = atomCoordinates[ii][jj]+getDeltaT()*velocities[ii][jj];
}
}
......@@ -188,7 +191,7 @@ void ReferenceStochasticDynamics::updatePart2( int numberOfAtoms, vector<RealVec
Update -- driver routine for performing stochastic dynamics update of coordinates
and velocities
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -196,7 +199,7 @@ void ReferenceStochasticDynamics::updatePart2( int numberOfAtoms, vector<RealVec
--------------------------------------------------------------------------------------- */
void ReferenceStochasticDynamics::update( int numberOfAtoms, vector<RealVec>& atomCoordinates,
void ReferenceStochasticDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses ){
// ---------------------------------------------------------------------------------------
......@@ -210,27 +213,15 @@ void ReferenceStochasticDynamics::update( int numberOfAtoms, vector<RealVec>& at
// first-time-through initialization
int numberOfAtoms = system.getNumParticles();
if( getTimeStep() == 0 ){
std::stringstream message;
message << methodName;
int errors = 0;
// invert masses
for( int ii = 0; ii < numberOfAtoms; ii++ ){
if( masses[ii] <= zero ){
message << "mass at atom index=" << ii << " (" << masses[ii] << ") is <= 0" << std::endl;
errors++;
} else {
inverseMasses[ii] = one/masses[ii];
}
}
// exit if errors
if( errors ){
SimTKOpenMMLog::printError( message );
if (masses[ii] == zero)
inverseMasses[ii] = zero;
else
inverseMasses[ii] = one/masses[ii];
}
}
......@@ -251,10 +242,12 @@ void ReferenceStochasticDynamics::update( int numberOfAtoms, vector<RealVec>& at
RealOpenMM invStepSize = 1.0/getDeltaT();
for (int i = 0; i < numberOfAtoms; ++i)
for (int j = 0; j < 3; ++j) {
velocities[i][j] = invStepSize*(xPrime[i][j]-atomCoordinates[i][j]);
atomCoordinates[i][j] = xPrime[i][j];
}
if (masses[i] != zero)
for (int j = 0; j < 3; ++j) {
velocities[i][j] = invStepSize*(xPrime[i][j]-atomCoordinates[i][j]);
atomCoordinates[i][j] = xPrime[i][j];
}
ReferenceVirtualSites::computePositions(system, atomCoordinates);
incrementTimeStep();
}
platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.h
View file @ fd01aa58
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -74,7 +74,7 @@ class ReferenceStochasticDynamics : public ReferenceDynamics {
Update
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -82,7 +82,7 @@ class ReferenceStochasticDynamics : public ReferenceDynamics {
--------------------------------------------------------------------------------------- */
void update( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses );
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/SimTKReference/ReferenceVariableStochasticDynamics.cpp
View file @ fd01aa58
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -29,6 +29,7 @@
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceVariableStochasticDynamics.h"
#include "ReferenceVirtualSites.h"
#include <cstdio>
......@@ -164,26 +165,13 @@ void ReferenceVariableStochasticDynamics::updatePart1( int numberOfAtoms, vector
// first-time-through initialization
if( getTimeStep() == 0 ){
std::stringstream message;
message << methodName;
int errors = 0;
// invert masses
for( int ii = 0; ii < numberOfAtoms; ii++ ){
if( masses[ii] <= 0 ){
message << "mass at atom index=" << ii << " (" << masses[ii] << ") is <= 0" << std::endl;
errors++;
} else {
inverseMasses[ii] = 1/masses[ii];
}
}
// exit if errors
if( errors ){
SimTKOpenMMLog::printError( message );
if (masses[ii] == 0)
inverseMasses[ii] = 0;
else
inverseMasses[ii] = 1/masses[ii];
}
}
......@@ -214,10 +202,12 @@ void ReferenceVariableStochasticDynamics::updatePart1( int numberOfAtoms, vector
const RealOpenMM noisescale = SQRT(2*kT/tau)*SQRT(0.5*(1-vscale*vscale)*tau);
for (int ii = 0; ii < numberOfAtoms; ii++) {
RealOpenMM sqrtInvMass = SQRT(inverseMasses[ii]);
for (int jj = 0; jj < 3; jj++) {
velocities[ii][jj] = vscale*velocities[ii][jj] + fscale*inverseMasses[ii]*forces[ii][jj] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
if (masses[ii] != 0) {
RealOpenMM sqrtInvMass = SQRT(inverseMasses[ii]);
for (int jj = 0; jj < 3; jj++) {
velocities[ii][jj] = vscale*velocities[ii][jj] + fscale*inverseMasses[ii]*forces[ii][jj] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
}
}
}
......@@ -248,9 +238,9 @@ void ReferenceVariableStochasticDynamics::updatePart2( int numberOfAtoms, vector
// perform second update
for (int ii = 0; ii < numberOfAtoms; ii++) {
for (int jj = 0; jj < 3; jj++) {
xPrime[ii][jj] = atomCoordinates[ii][jj]+getDeltaT()*velocities[ii][jj];
}
if (inverseMasses[ii] != 0.0)
for (int jj = 0; jj < 3; jj++)
xPrime[ii][jj] = atomCoordinates[ii][jj]+getDeltaT()*velocities[ii][jj];
}
}
......@@ -260,7 +250,7 @@ void ReferenceVariableStochasticDynamics::updatePart2( int numberOfAtoms, vector
Update -- driver routine for performing stochastic dynamics update of coordinates
and velocities
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -268,7 +258,7 @@ void ReferenceVariableStochasticDynamics::updatePart2( int numberOfAtoms, vector
--------------------------------------------------------------------------------------- */
void ReferenceVariableStochasticDynamics::update( int numberOfAtoms, vector<RealVec>& atomCoordinates,
void ReferenceVariableStochasticDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM maxStepSize ){
......@@ -280,6 +270,7 @@ void ReferenceVariableStochasticDynamics::update( int numberOfAtoms, vector<Real
// 1st update
int numberOfAtoms = system.getNumParticles();
updatePart1( numberOfAtoms, atomCoordinates, velocities, forces, masses, inverseMasses, xPrime, maxStepSize );
// 2nd update
......@@ -294,12 +285,14 @@ void ReferenceVariableStochasticDynamics::update( int numberOfAtoms, vector<Real
// copy xPrime -> atomCoordinates
for( int ii = 0; ii < numberOfAtoms; ii++ ){
atomCoordinates[ii][0] = xPrime[ii][0];
atomCoordinates[ii][1] = xPrime[ii][1];
atomCoordinates[ii][2] = xPrime[ii][2];
for( int ii = 0; ii < numberOfAtoms; ii++ ) {
if (masses[ii] != 0.0) {
atomCoordinates[ii][0] = xPrime[ii][0];
atomCoordinates[ii][1] = xPrime[ii][1];
atomCoordinates[ii][2] = xPrime[ii][2];
}
}
ReferenceVirtualSites::computePositions(system, atomCoordinates);
incrementTimeStep();
}
platforms/reference/src/SimTKReference/ReferenceVariableStochasticDynamics.h
View file @ fd01aa58
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -92,7 +92,7 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
Update
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -101,7 +101,7 @@ class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
--------------------------------------------------------------------------------------- */
void update( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses, RealOpenMM maxStepSize );
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/SimTKReference/ReferenceVariableVerletDynamics.cpp
View file @ fd01aa58
/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Peter Eastman, Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -29,6 +29,7 @@
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceVariableVerletDynamics.h"
#include "ReferenceVirtualSites.h"
using std::vector;
using OpenMM::RealVec;
......@@ -45,18 +46,8 @@ using OpenMM::RealVec;
ReferenceVariableVerletDynamics::ReferenceVariableVerletDynamics( int numberOfAtoms, RealOpenMM accuracy ) :
ReferenceDynamics( numberOfAtoms, 0.0f, 0.0f ), _accuracy(accuracy) {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceVariableVerletDynamics::ReferenceVariableVerletDynamics";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
xPrime.resize(numberOfAtoms);
inverseMasses.resize(numberOfAtoms);
xPrime.resize(numberOfAtoms);
inverseMasses.resize(numberOfAtoms);
}
/**---------------------------------------------------------------------------------------
......@@ -102,7 +93,7 @@ void ReferenceVariableVerletDynamics::setAccuracy( RealOpenMM accuracy ) {
Update -- driver routine for performing Verlet dynamics update of coordinates
and velocities
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -111,7 +102,7 @@ void ReferenceVariableVerletDynamics::setAccuracy( RealOpenMM accuracy ) {
--------------------------------------------------------------------------------------- */
void ReferenceVariableVerletDynamics::update( int numberOfAtoms, vector<RealVec>& atomCoordinates,
void ReferenceVariableVerletDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM maxStepSize ){
......@@ -126,27 +117,15 @@ void ReferenceVariableVerletDynamics::update( int numberOfAtoms, vector<RealVec>
// first-time-through initialization
int numberOfAtoms = system.getNumParticles();
if( getTimeStep() == 0 ){
std::stringstream message;
message << methodName;
int errors = 0;
// invert masses
for( int ii = 0; ii < numberOfAtoms; ii++ ){
if( masses[ii] <= zero ){
message << "mass at atom index=" << ii << " (" << masses[ii] << ") is <= 0" << std::endl;
errors++;
} else {
inverseMasses[ii] = one/masses[ii];
}
}
// exit if errors
if( errors ){
SimTKOpenMMLog::printError( message );
if (masses[ii] == zero)
inverseMasses[ii] = zero;
else
inverseMasses[ii] = one/masses[ii];
}
}
......@@ -168,10 +147,11 @@ void ReferenceVariableVerletDynamics::update( int numberOfAtoms, vector<RealVec>
RealOpenMM vstep = 0.5f*(newStepSize+getDeltaT()); // The time interval by which to advance the velocities
setDeltaT(newStepSize);
for (int i = 0; i < numberOfAtoms; ++i) {
for (int j = 0; j < 3; ++j) {
RealOpenMM vPrime = velocities[i][j] + inverseMasses[i]*forces[i][j]*vstep;
xPrime[i][j] = atomCoordinates[i][j] + vPrime*getDeltaT();
}
if (masses[i] != zero)
for (int j = 0; j < 3; ++j) {
RealOpenMM vPrime = velocities[i][j] + inverseMasses[i]*forces[i][j]*vstep;
xPrime[i][j] = atomCoordinates[i][j] + vPrime*getDeltaT();
}
}
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if (referenceConstraintAlgorithm)
......@@ -181,12 +161,13 @@ void ReferenceVariableVerletDynamics::update( int numberOfAtoms, vector<RealVec>
RealOpenMM velocityScale = one/getDeltaT();
for (int i = 0; i < numberOfAtoms; ++i) {
for (int j = 0; j < 3; ++j) {
velocities[i][j] = velocityScale*(xPrime[i][j] - atomCoordinates[i][j]);
atomCoordinates[i][j] = xPrime[i][j];
}
if (masses[i] != zero)
for (int j = 0; j < 3; ++j) {
velocities[i][j] = velocityScale*(xPrime[i][j] - atomCoordinates[i][j]);
atomCoordinates[i][j] = xPrime[i][j];
}
}
ReferenceVirtualSites::computePositions(system, atomCoordinates);
incrementTimeStep();
}
......
platforms/reference/src/SimTKReference/ReferenceVariableVerletDynamics.h
View file @ fd01aa58
/* Portions copyright (c) 2006-2009 Stanford University and Simbios.
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Peter Eastman, Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -80,7 +80,7 @@ class ReferenceVariableVerletDynamics : public ReferenceDynamics {
Update
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -89,7 +89,7 @@ class ReferenceVariableVerletDynamics : public ReferenceDynamics {
--------------------------------------------------------------------------------------- */
void update( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses, RealOpenMM maxStepSize );
};
......
platforms/reference/src/SimTKReference/ReferenceVerletDynamics.cpp
View file @ fd01aa58
/* Portions copyright (c) 2006-2008 Stanford University and Simbios.
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Peter Eastman, Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -29,6 +29,7 @@
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceVerletDynamics.h"
#include "ReferenceVirtualSites.h"
#include <cstdio>
......@@ -84,7 +85,7 @@ ReferenceVerletDynamics::~ReferenceVerletDynamics( ){
Update -- driver routine for performing Verlet dynamics update of coordinates
and velocities
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -92,7 +93,7 @@ ReferenceVerletDynamics::~ReferenceVerletDynamics( ){
--------------------------------------------------------------------------------------- */
void ReferenceVerletDynamics::update( int numberOfAtoms, vector<RealVec>& atomCoordinates,
void ReferenceVerletDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& masses ){
......@@ -107,37 +108,26 @@ void ReferenceVerletDynamics::update( int numberOfAtoms, vector<RealVec>& atomCo
// first-time-through initialization
int numberOfAtoms = system.getNumParticles();
if( getTimeStep() == 0 ){
std::stringstream message;
message << methodName;
int errors = 0;
// invert masses
for( int ii = 0; ii < numberOfAtoms; ii++ ){
if( masses[ii] <= zero ){
message << "mass at atom index=" << ii << " (" << masses[ii] << ") is <= 0" << std::endl;
errors++;
} else {
inverseMasses[ii] = one/masses[ii];
}
}
// exit if errors
if( errors ){
SimTKOpenMMLog::printError( message );
if (masses[ii] == zero)
inverseMasses[ii] = zero;
else
inverseMasses[ii] = one/masses[ii];
}
}
// Perform the integration.
for (int i = 0; i < numberOfAtoms; ++i) {
for (int j = 0; j < 3; ++j) {
velocities[i][j] += inverseMasses[i]*forces[i][j]*getDeltaT();
xPrime[i][j] = atomCoordinates[i][j] + velocities[i][j]*getDeltaT();
}
if (masses[i] != zero)
for (int j = 0; j < 3; ++j) {
velocities[i][j] += inverseMasses[i]*forces[i][j]*getDeltaT();
xPrime[i][j] = atomCoordinates[i][j] + velocities[i][j]*getDeltaT();
}
}
ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
if( referenceConstraintAlgorithm )
......@@ -147,11 +137,13 @@ void ReferenceVerletDynamics::update( int numberOfAtoms, vector<RealVec>& atomCo
RealOpenMM velocityScale = static_cast<RealOpenMM>( 1.0/getDeltaT() );
for (int i = 0; i < numberOfAtoms; ++i) {
for (int j = 0; j < 3; ++j) {
velocities[i][j] = velocityScale*(xPrime[i][j] - atomCoordinates[i][j]);
atomCoordinates[i][j] = xPrime[i][j];
}
if (masses[i] != zero)
for (int j = 0; j < 3; ++j) {
velocities[i][j] = velocityScale*(xPrime[i][j] - atomCoordinates[i][j]);
atomCoordinates[i][j] = xPrime[i][j];
}
}
ReferenceVirtualSites::computePositions(system, atomCoordinates);
incrementTimeStep();
}
platforms/reference/src/SimTKReference/ReferenceVerletDynamics.h
View file @ fd01aa58
/* Portions copyright (c) 2006-2008 Stanford University and Simbios.
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -63,7 +63,7 @@ class ReferenceVerletDynamics : public ReferenceDynamics {
Update
@param numberOfAtoms number of atoms
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
......@@ -71,7 +71,7 @@ class ReferenceVerletDynamics : public ReferenceDynamics {
--------------------------------------------------------------------------------------- */
void update( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
void update(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses );
};
......
platforms/reference/src/SimTKReference/ReferenceVirtualSites.cpp 0 → 100644
View file @ fd01aa58
/* -------------------------------------------------------------------------- *
* 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) 2012 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 "ReferenceVirtualSites.h"
#include "openmm/VirtualSite.h"
using namespace OpenMM;
using namespace std;
void ReferenceVirtualSites::computePositions(const OpenMM::System& system, vector<OpenMM::RealVec>& atomCoordinates) {
for (int i = 0; i < system.getNumParticles(); i++)
if (system.isVirtualSite(i)) {
if (dynamic_cast<const VirtualSite::TwoParticleAverage*>(&system.getVirtualSite(i)) != NULL) {
// A two particle average.
const VirtualSite::TwoParticleAverage& site = dynamic_cast<const VirtualSite::TwoParticleAverage&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1);
RealOpenMM w1 = site.getWeight(0), w2 = site.getWeight(1);
atomCoordinates[i] = atomCoordinates[p1]*w1 + atomCoordinates[p2]*w2;
}
else if (dynamic_cast<const VirtualSite::ThreeParticleAverage*>(&system.getVirtualSite(i)) != NULL) {
// A three particle average.
const VirtualSite::ThreeParticleAverage& site = dynamic_cast<const VirtualSite::ThreeParticleAverage&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
RealOpenMM w1 = site.getWeight(0), w2 = site.getWeight(1), w3 = site.getWeight(2);
atomCoordinates[i] = atomCoordinates[p1]*w1 + atomCoordinates[p2]*w2 + atomCoordinates[p3]*w3;
}
else if (dynamic_cast<const VirtualSite::OutOfPlane*>(&system.getVirtualSite(i)) != NULL) {
// An out of plane site.
const VirtualSite::OutOfPlane& site = dynamic_cast<const VirtualSite::OutOfPlane&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
RealOpenMM w12 = site.getWeight12(), w13 = site.getWeight13(), wcross = site.getWeightCross();
RealVec v12 = atomCoordinates[p2]-atomCoordinates[p1];
RealVec v13 = atomCoordinates[p3]-atomCoordinates[p1];
RealVec cross = v12.cross(v13);
atomCoordinates[i] = atomCoordinates[p1] + v12*w12 + v13*w13 + cross*wcross;
}
}
}
void ReferenceVirtualSites::distributeForces(const OpenMM::System& system, const vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& forces) {
for (int i = 0; i < system.getNumParticles(); i++)
if (system.isVirtualSite(i)) {
RealVec f = forces[i];
if (dynamic_cast<const VirtualSite::TwoParticleAverage*>(&system.getVirtualSite(i)) != NULL) {
// A two particle average.
const VirtualSite::TwoParticleAverage& site = dynamic_cast<const VirtualSite::TwoParticleAverage&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1);
RealOpenMM w1 = site.getWeight(0), w2 = site.getWeight(1);
forces[p1] += f*w1;
forces[p2] += f*w2;
}
else if (dynamic_cast<const VirtualSite::ThreeParticleAverage*>(&system.getVirtualSite(i)) != NULL) {
// A three particle average.
const VirtualSite::ThreeParticleAverage& site = dynamic_cast<const VirtualSite::ThreeParticleAverage&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
RealOpenMM w1 = site.getWeight(0), w2 = site.getWeight(1), w3 = site.getWeight(2);
forces[p1] += f*w1;
forces[p2] += f*w2;
forces[p3] += f*w3;
}
else if (dynamic_cast<const VirtualSite::OutOfPlane*>(&system.getVirtualSite(i)) != NULL) {
// An out of plane site.
const VirtualSite::OutOfPlane& site = dynamic_cast<const VirtualSite::OutOfPlane&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
RealOpenMM w12 = site.getWeight12(), w13 = site.getWeight13(), wcross = site.getWeightCross();
RealVec v12 = atomCoordinates[p2]-atomCoordinates[p1];
RealVec v13 = atomCoordinates[p3]-atomCoordinates[p1];
RealVec f2(w12*f[0] - wcross*v13[2]*f[1] + wcross*v13[1]*f[2],
wcross*v13[2]*f[0] + w12*f[1] - wcross*v13[0]*f[2],
-wcross*v13[1]*f[0] + wcross*v13[0]*f[1] + w12*f[2]);
RealVec f3(w13*f[0] + wcross*v12[2]*f[1] - wcross*v12[1]*f[2],
-wcross*v12[2]*f[0] + w13*f[1] + wcross*v12[0]*f[2],
wcross*v12[1]*f[0] - wcross*v12[0]*f[1] + w13*f[2]);
forces[p1] += f-f2-f3;
forces[p2] += f2;
forces[p3] += f3;
}
}
}
platforms/reference/src/SimTKReference/ReferenceVirtualSites.h 0 → 100644
View file @ fd01aa58
/* -------------------------------------------------------------------------- *
* 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) 2012 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. *
* -------------------------------------------------------------------------- */
#ifndef __ReferenceVirtualSites_H__
#define __ReferenceVirtualSites_H__
#include "openmm/System.h"
#include "../SimTKUtilities/RealVec.h"
#include <vector>
class ReferenceVirtualSites {
private:
const OpenMM::System& system;
public:
/**
* Compute the positions of all virtual sites.
*/
static void computePositions(const OpenMM::System& system, std::vector<OpenMM::RealVec>& atomCoordinates);
/**
* Distribute forces from virtual sites to the atoms they are based on.
*/
static void distributeForces(const OpenMM::System& system, const std::vector<OpenMM::RealVec>& atomCoordinates, std::vector<OpenMM::RealVec>& forces);
};
#endif // __ReferenceVirtualSites_H__
platforms/reference/tests/TestReferenceVirtualSites.cpp 0 → 100644
View file @ fd01aa58
/* -------------------------------------------------------------------------- *
* 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) 2012 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 virtual sites.
*/
#include "../../../tests/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/CustomBondForce.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/VirtualSite.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
/**
* Check that massless particles are handled correctly.
*/
void testMasslessParticle() {
ReferencePlatform platform;
System system;
system.addParticle(0.0);
system.addParticle(1.0);
CustomBondForce* bonds = new CustomBondForce("-1/r");
system.addForce(bonds);
vector<double> params;
bonds->addBond(0, 1, params);
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
vector<Vec3> velocities(2);
velocities[0] = Vec3(0, 0, 0);
velocities[1] = Vec3(0, 1, 0);
context.setVelocities(velocities);
// The second particle should move in a circular orbit around the first one.
// Compare it to the analytical solution.
for (int i = 0; i < 1000; ++i) {
State state = context.getState(State::Positions | State::Velocities | State::Forces);
double time = state.getTime();
ASSERT_EQUAL_VEC(Vec3(), state.getPositions()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(), state.getVelocities()[0], 0.0);
ASSERT_EQUAL_VEC(Vec3(cos(time), sin(time), 0), state.getPositions()[1], 0.01);
ASSERT_EQUAL_VEC(Vec3(-sin(time), cos(time), 0), state.getVelocities()[1], 0.01);
integrator.step(1);
}
}
/**
* Test a TwoParticleAverage virtual site.
*/
void testTwoParticleAverage() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new VirtualSite::TwoParticleAverage(0, 1, 0.8, 0.2));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.8+pos[1]*0.2, pos[2], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.1+0.3*0.8, 0, 0), state.getForces()[0], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.2+0.3*0.2, 0, 0), state.getForces()[1], 1e-10);
}
}
/**
* Test a ThreeParticleAverage virtual site.
*/
void testThreeParticleAverage() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new VirtualSite::ThreeParticleAverage(0, 1, 2, 0.2, 0.3, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
ASSERT_EQUAL_VEC(pos[0]*0.2+pos[1]*0.3+pos[2]*0.5, pos[3], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4*0.2, 0, 0), state.getForces()[0], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.2+0.4*0.3, 0, 0), state.getForces()[1], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.3+0.4*0.5, 0, 0), state.getForces()[2], 1e-10);
}
}
/**
* Test an OutOfPlane virtual site.
*/
void testOutOfPlane() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new VirtualSite::OutOfPlane(0, 1, 2, 0.3, 0.4, 0.5));
CustomExternalForce* forceField = new CustomExternalForce("-a*x");
system.addForce(forceField);
forceField->addPerParticleParameter("a");
vector<double> params(1);
params[0] = 0.1;
forceField->addParticle(0, params);
params[0] = 0.2;
forceField->addParticle(1, params);
params[0] = 0.3;
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
context.setPositions(positions);
context.applyConstraints(0.0001);
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 v12 = pos[1]-pos[0];
Vec3 v13 = pos[2]-pos[0];
Vec3 cross = v12.cross(v13);
ASSERT_EQUAL_VEC(pos[0]+v12*0.3+v13*0.4+cross*0.5, pos[3], 1e-10);
const vector<Vec3>& f = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0.1+0.2+0.3+0.4, 0, 0), f[0]+f[1]+f[2], 1e-10);
Vec3 f2(0.4*0.3, 0.4*0.5*v13[2], -0.4*0.5*v13[1]);
Vec3 f3(0.4*0.4, -0.4*0.5*v12[2], 0.4*0.5*v12[1]);
ASSERT_EQUAL_VEC(Vec3(0.1+0.4, 0, 0)-f2-f3, f[0], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.2, 0, 0)+f2, f[1], 1e-10);
ASSERT_EQUAL_VEC(Vec3(0.3, 0, 0)+f3, f[2], 1e-10);
}
}
/**
* Make sure that energy, linear momentum, and angular momentum are all conserved
* when using virtual sites.
*/
void testConservationLaws() {
ReferencePlatform platform;
System system;
NonbondedForce* forceField = new NonbondedForce();
system.addForce(forceField);
vector<Vec3> positions;
// Create a linear molecule with a TwoParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(2, new VirtualSite::TwoParticleAverage(0, 1, 0.4, 0.6));
system.addConstraint(0, 1, 2.0);
for (int i = 0; i < 3; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i, j, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(2, 0, 0));
positions.push_back(Vec3());
// Create a planar molecule with a ThreeParticleAverage virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(6, new VirtualSite::ThreeParticleAverage(3, 4, 5, 0.3, 0.5, 0.2));
system.addConstraint(3, 4, 1.0);
system.addConstraint(3, 5, 1.0);
system.addConstraint(4, 5, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+3, j+3, 0, 1, 0);
}
positions.push_back(Vec3(0, 0, 1));
positions.push_back(Vec3(1, 0, 1));
positions.push_back(Vec3(0, 1, 1));
positions.push_back(Vec3());
// Create a tetrahedral molecule with an OutOfPlane virtual site.
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(10, new VirtualSite::OutOfPlane(7, 8, 9, 0.3, 0.5, 0.2));
system.addConstraint(7, 8, 1.0);
system.addConstraint(7, 9, 1.0);
system.addConstraint(8, 9, sqrt(2.0));
for (int i = 0; i < 4; i++) {
forceField->addParticle(0, 1, 10);
for (int j = 0; j < i; j++)
forceField->addException(i+7, j+7, 0, 1, 0);
}
positions.push_back(Vec3(1, 0, -1));
positions.push_back(Vec3(2, 0, -1));
positions.push_back(Vec3(1, 1, -1));
positions.push_back(Vec3());
// Simulate it and check conservation laws.
VerletIntegrator integrator(0.002);
Context context(system, integrator, platform);
context.setPositions(positions);
context.applyConstraints(0.0001);
int numParticles = system.getNumParticles();
double initialEnergy;
Vec3 initialMomentum, initialAngularMomentum;
for (int i = 0; i < 1000; i++) {
State state = context.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
const vector<Vec3>& pos = state.getPositions();
const vector<Vec3>& vel = state.getVelocities();
const vector<Vec3>& f = state.getForces();
double energy = state.getPotentialEnergy();
for (int i = 0; i < numParticles; i++) {
Vec3 v = vel[i] + f[i]*0.5*integrator.getStepSize();
energy += 0.5*system.getParticleMass(i)*v.dot(v);
}
if (i == 0)
initialEnergy = energy;
else
ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
Vec3 momentum;
for (int j = 0; j < numParticles; j++)
momentum += vel[j]*system.getParticleMass(j);
if (i == 0)
initialMomentum = momentum;
else
ASSERT_EQUAL_VEC(initialMomentum, momentum, 1e-10);
Vec3 angularMomentum;
for (int j = 0; j < numParticles; j++)
angularMomentum += pos[j].cross(vel[j])*system.getParticleMass(j);
if (i == 0)
initialAngularMomentum = angularMomentum;
else
ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 1e-10);
integrator.step(1);
}
}
int main() {
try {
testMasslessParticle();
testTwoParticleAverage();
testThreeParticleAverage();
testOutOfPlane();
testConservationLaws();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
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