Skip to content

GitLab

Commit 94fbbe9c authored by Peter Eastman's avatar Peter Eastman
Browse files

Reference implementation of periodic boundary conditions for AMOEBA bonded forces

parent b4dcef47
Hide whitespace changes
Inline Side-by-side
Showing with 549 additions and 209 deletions
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceBondForce.cpp
View file @ 94fbbe9c
......@@ -28,6 +28,13 @@
using std::vector;
using namespace OpenMM;
void AmoebaReferenceBondForce::setPeriodic(OpenMM::RealVec* vectors) {
usePeriodic = true;
boxVectors[0] = vectors[0];
boxVectors[1] = vectors[1];
boxVectors[2] = vectors[2];
}
/**---------------------------------------------------------------------------------------
Calculate Amoeba bond ixn (force and energy)
......@@ -63,7 +70,10 @@ RealOpenMM AmoebaReferenceBondForce::calculateBondIxn(const RealVec& positionAto
// get deltaR, R2, and R between 2 atoms
std::vector<RealOpenMM> deltaR;
AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomB, deltaR);
if (usePeriodic)
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomA, positionAtomB, deltaR, boxVectors);
else
AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomB, deltaR);
RealOpenMM r = AmoebaReferenceForce::getNorm3(deltaR);
// deltaIdeal = r - r_0
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceBondForce.h
View file @ 94fbbe9c
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -40,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceBondForce() {};
AmoebaReferenceBondForce() : usePeriodic(false) {};
/**---------------------------------------------------------------------------------------
......@@ -49,7 +49,16 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferenceBondForce() {};
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions.
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
......@@ -79,6 +88,9 @@ public:
private:
bool usePeriodic;
RealVec boxVectors[3];
/**---------------------------------------------------------------------------------------
Calculate Amoeba bond ixns (force and energy)
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceForce.cpp
View file @ 94fbbe9c
......@@ -39,18 +39,33 @@ using namespace OpenMM;
void AmoebaReferenceForce::loadDeltaR(const RealVec& xVector, const RealVec& yVector,
std::vector<RealOpenMM>& deltaR) {
deltaR.resize(0);
deltaR.push_back(yVector[0] - xVector[0]);
deltaR.push_back(yVector[1] - xVector[1]);
deltaR.push_back(yVector[2] - xVector[2]);
}
/**---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
Load delta of two vectors, applying periodic boundary conditions
//static const std::string methodName = "AmoebaReferenceForce::loadDeltaR";
@param xVector first vector
@param yVector second vector
@param deltaR output vector: y - x
@param boxVectors periodic box vectors
// ---------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------- */
void AmoebaReferenceForce::loadDeltaRPeriodic(const RealVec& xVector, const RealVec& yVector,
std::vector<RealOpenMM>& deltaR, const RealVec* boxVectors) {
RealVec diff = yVector-xVector;
diff -= boxVectors[2]*floor(diff[2]/boxVectors[2][2]+0.5);
diff -= boxVectors[1]*floor(diff[1]/boxVectors[1][1]+0.5);
diff -= boxVectors[0]*floor(diff[0]/boxVectors[0][0]+0.5);
deltaR.resize(0);
deltaR.push_back(yVector[0] - xVector[0]);
deltaR.push_back(yVector[1] - xVector[1]);
deltaR.push_back(yVector[2] - xVector[2]);
deltaR.push_back(diff[0]);
deltaR.push_back(diff[1]);
deltaR.push_back(diff[2]);
}
/**---------------------------------------------------------------------------------------
......@@ -64,13 +79,6 @@ void AmoebaReferenceForce::loadDeltaR(const RealVec& xVector, const RealVec& yVe
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceForce::getNormSquared3(const std::vector<RealOpenMM>& inputVector) {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getNorm3";
// ---------------------------------------------------------------------------------------
// get 3 norm
return (inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2]);
......@@ -87,13 +95,6 @@ RealOpenMM AmoebaReferenceForce::getNormSquared3(const std::vector<RealOpenMM>&
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceForce::getNormSquared3(const RealOpenMM* inputVector) {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getNorm3";
// ---------------------------------------------------------------------------------------
// get 3 norm
return (inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2]);
......@@ -110,39 +111,18 @@ RealOpenMM AmoebaReferenceForce::getNormSquared3(const RealOpenMM* inputVector)
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceForce::getNorm3(const std::vector<RealOpenMM>& inputVector) {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getNorm3";
// ---------------------------------------------------------------------------------------
// get 3 norm
return SQRT(inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2]);
}
RealOpenMM AmoebaReferenceForce::getNorm3(const RealOpenMM* inputVector) {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getNorm3";
// ---------------------------------------------------------------------------------------
// get 3 norm
return SQRT(inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2]);
}
RealOpenMM AmoebaReferenceForce::normalizeVector3(RealOpenMM* inputVector) {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::normalizeVector3";
// ---------------------------------------------------------------------------------------
RealOpenMM norm = SQRT(inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2]);
if (norm > 0.0) {
RealOpenMM normI = 1.0/norm;
......@@ -166,13 +146,6 @@ RealOpenMM AmoebaReferenceForce::normalizeVector3(RealOpenMM* inputVector) {
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceForce::getDotProduct3(const std::vector<RealOpenMM>& xVector, const std::vector<RealOpenMM>& yVector) {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getDotProduct3";
// ---------------------------------------------------------------------------------------
// get dot product
return xVector[0]*yVector[0] + xVector[1]*yVector[1] + xVector[2]*yVector[2];
......@@ -190,26 +163,12 @@ RealOpenMM AmoebaReferenceForce::getDotProduct3(const std::vector<RealOpenMM>& x
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceForce::getDotProduct3(const RealOpenMM* xVector, const RealOpenMM* yVector) {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getDotProduct3";
// ---------------------------------------------------------------------------------------
// get dot product
return xVector[0]*yVector[0] + xVector[1]*yVector[1] + xVector[2]*yVector[2];
}
RealOpenMM AmoebaReferenceForce::getDotProduct3(const RealOpenMM* xVector, const OpenMM::Vec3& yVector) {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getDotProduct3";
// ---------------------------------------------------------------------------------------
// get dot product
return xVector[0]*yVector[0] + xVector[1]*yVector[1] + xVector[2]*yVector[2];
......@@ -228,13 +187,6 @@ RealOpenMM AmoebaReferenceForce::getDotProduct3(const RealOpenMM* xVector, const
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceForce::getDotProduct3(unsigned int vectorOffset, const std::vector<RealOpenMM>& xVector, const RealOpenMM* yVector) {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getDotProduct3";
// ---------------------------------------------------------------------------------------
// get dot product
return xVector[vectorOffset+0]*yVector[0] + xVector[vectorOffset+1]*yVector[1] + xVector[vectorOffset+2]*yVector[2];
......@@ -253,18 +205,9 @@ RealOpenMM AmoebaReferenceForce::getDotProduct3(unsigned int vectorOffset, const
void AmoebaReferenceForce::getCrossProduct(const std::vector<RealOpenMM>& xVector,
const std::vector<RealOpenMM>& yVector,
std::vector<RealOpenMM>& zVector) {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getCrossProduct";
// ---------------------------------------------------------------------------------------
zVector[0] = xVector[1]*yVector[2] - xVector[2]*yVector[1];
zVector[1] = xVector[2]*yVector[0] - xVector[0]*yVector[2];
zVector[2] = xVector[0]*yVector[1] - xVector[1]*yVector[0];
return;
}
/**---------------------------------------------------------------------------------------
......@@ -280,17 +223,8 @@ void AmoebaReferenceForce::getCrossProduct(const std::vector<RealOpenMM>& xVecto
void AmoebaReferenceForce::getCrossProduct(const RealOpenMM* xVector,
const RealOpenMM* yVector,
RealOpenMM* zVector) {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceForce::getCrossProduct";
// ---------------------------------------------------------------------------------------
zVector[0] = xVector[1]*yVector[2] - xVector[2]*yVector[1];
zVector[1] = xVector[2]*yVector[0] - xVector[0]*yVector[2];
zVector[2] = xVector[0]*yVector[1] - xVector[1]*yVector[0];
return;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceForce.h
View file @ 94fbbe9c
......@@ -64,6 +64,18 @@ public:
static void loadDeltaR(const OpenMM::RealVec& xVector, const OpenMM::RealVec& yVector,
std::vector<RealOpenMM>& deltaR);
/**---------------------------------------------------------------------------------------
Load delta of two vectors, applying periodic boundary conditions
@param xVector first vector
@param yVector second vector
@param deltaR output vector: y - x
@param boxVectors periodic box vectors
--------------------------------------------------------------------------------------- */
static void loadDeltaRPeriodic(const RealVec& xVector, const RealVec& yVector, std::vector<RealOpenMM>& deltaR, const RealVec* boxVectors);
/**---------------------------------------------------------------------------------------
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceInPlaneAngleForce.cpp
View file @ 94fbbe9c
......@@ -24,10 +24,18 @@
#include "AmoebaReferenceForce.h"
#include "AmoebaReferenceInPlaneAngleForce.h"
#include "ReferenceForce.h"
using std::vector;
using namespace OpenMM;
void AmoebaReferenceInPlaneAngleForce::setPeriodic(OpenMM::RealVec* vectors) {
usePeriodic = true;
boxVectors[0] = vectors[0];
boxVectors[1] = vectors[1];
boxVectors[2] = vectors[2];
}
/**---------------------------------------------------------------------------------------
Get dEdT and energy prefactor given cosine of angle :: the calculation for different
......@@ -120,17 +128,6 @@ RealOpenMM AmoebaReferenceInPlaneAngleForce::calculateAngleIxn(const RealVec& po
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealVec* forces) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceInPlaneAngleForce::calculateAngleIxn";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
// ---------------------------------------------------------------------------------------
// T = AD x CD
// P = B + T*delta
// AP = A - P
......@@ -139,40 +136,48 @@ RealOpenMM AmoebaReferenceInPlaneAngleForce::calculateAngleIxn(const RealVec& po
enum { AD, BD, CD, T, AP, P, CP, M, APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex };
std::vector<RealOpenMM> deltaR[LastDeltaAtomIndex];
for (int ii = 0; ii < LastDeltaAtomIndex; ii++) {
deltaR[ii].resize(3);
RealVec deltaR[LastDeltaAtomIndex];
if (usePeriodic) {
deltaR[AD] = ReferenceForce::getDeltaRPeriodic(positionAtomD, positionAtomA, boxVectors);
deltaR[BD] = ReferenceForce::getDeltaRPeriodic(positionAtomD, positionAtomB, boxVectors);
deltaR[CD] = ReferenceForce::getDeltaRPeriodic(positionAtomD, positionAtomC, boxVectors);
}
else {
deltaR[AD] = ReferenceForce::getDeltaR(positionAtomD, positionAtomA);
deltaR[BD] = ReferenceForce::getDeltaR(positionAtomD, positionAtomB);
deltaR[CD] = ReferenceForce::getDeltaR(positionAtomD, positionAtomC);
}
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomA, deltaR[AD]);
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomB, deltaR[BD]);
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomC, deltaR[CD]);
AmoebaReferenceForce::getCrossProduct(deltaR[AD], deltaR[CD], deltaR[T]);
deltaR[T] = deltaR[AD].cross(deltaR[CD]);
RealOpenMM rT2 = AmoebaReferenceForce::getNormSquared3(deltaR[T]);
RealOpenMM delta = AmoebaReferenceForce::getDotProduct3(deltaR[T], deltaR[BD])/rT2;
delta *= -one;
for (int ii = 0; ii < 3; ii++) {
deltaR[P][ii] = positionAtomB[ii] + deltaR[T][ii]*delta;
deltaR[AP][ii] = positionAtomA[ii] - deltaR[P][ii];
deltaR[CP][ii] = positionAtomC[ii] - deltaR[P][ii];
RealOpenMM rT2 = deltaR[T].dot(deltaR[T]);
RealOpenMM delta = deltaR[T].dot(deltaR[BD])/rT2;
delta *= -1;
deltaR[P] = positionAtomB + deltaR[T]*delta;
if (usePeriodic) {
deltaR[AP] = ReferenceForce::getDeltaRPeriodic(deltaR[P], positionAtomA, boxVectors);
deltaR[CP] = ReferenceForce::getDeltaRPeriodic(deltaR[P], positionAtomC, boxVectors);
}
else {
deltaR[AP] = ReferenceForce::getDeltaR(deltaR[P], positionAtomA);
deltaR[CP] = ReferenceForce::getDeltaR(deltaR[P], positionAtomC);
}
RealOpenMM rAp2 = AmoebaReferenceForce::getNormSquared3(deltaR[AP]);
RealOpenMM rCp2 = AmoebaReferenceForce::getNormSquared3(deltaR[CP]);
if (rAp2 <= zero && rCp2 <= zero) {
return zero;
RealOpenMM rAp2 = deltaR[AP].dot(deltaR[AP]);
RealOpenMM rCp2 = deltaR[CP].dot(deltaR[CP]);
if (rAp2 <= 0 && rCp2 <= 0) {
return 0;
}
AmoebaReferenceForce::getCrossProduct(deltaR[CP], deltaR[AP], deltaR[M]);
deltaR[M] = deltaR[CP].cross(deltaR[AP]);
RealOpenMM rm = AmoebaReferenceForce::getNorm3(deltaR[M]);
RealOpenMM rm = SQRT(deltaR[M].dot(deltaR[M]));
if (rm < 1.0e-06) {
rm = 1.0e-06;
}
RealOpenMM dot = AmoebaReferenceForce::getDotProduct3(deltaR[AP], deltaR[CP]);
RealOpenMM dot = deltaR[AP].dot(deltaR[CP]);
RealOpenMM cosine = dot/SQRT(rAp2*rCp2);
RealOpenMM dEdR;
......@@ -182,64 +187,47 @@ RealOpenMM AmoebaReferenceInPlaneAngleForce::calculateAngleIxn(const RealVec& po
RealOpenMM termA = -dEdR/(rAp2*rm);
RealOpenMM termC = dEdR/(rCp2*rm);
AmoebaReferenceForce::getCrossProduct(deltaR[AP], deltaR[M], deltaR[APxM]);
AmoebaReferenceForce::getCrossProduct(deltaR[CP], deltaR[M], deltaR[CPxM]);
deltaR[APxM] = deltaR[AP].cross(deltaR[M]);
deltaR[CPxM] = deltaR[CP].cross(deltaR[M]);
// forces will be gathered here
enum { dA, dB, dC, dD, LastDIndex };
std::vector<RealOpenMM> forceTerm[LastDIndex];
for (int ii = 0; ii < LastDIndex; ii++) {
forceTerm[ii].resize(3);
}
RealVec forceTerm[LastDIndex];
for (int ii = 0; ii < 3; ii++) {
forceTerm[dA][ii] = deltaR[APxM][ii]*termA;
forceTerm[dC][ii] = deltaR[CPxM][ii]*termC;
forceTerm[dB][ii] = -one*(forceTerm[dA][ii] + forceTerm[dC][ii]);
}
RealOpenMM pTrT2 = AmoebaReferenceForce::getDotProduct3(forceTerm[dB], deltaR[T]);
forceTerm[dA] = deltaR[APxM]*termA;
forceTerm[dC] = deltaR[CPxM]*termC;
forceTerm[dB] = -(forceTerm[dA] + forceTerm[dC]);
RealOpenMM pTrT2 = forceTerm[dB].dot(deltaR[T]);
pTrT2 /= rT2;
AmoebaReferenceForce::getCrossProduct(deltaR[CD], forceTerm[dB], deltaR[CDxdB]);
AmoebaReferenceForce::getCrossProduct(forceTerm[dB], deltaR[AD], deltaR[dBxAD]);
deltaR[CDxdB] = deltaR[CD].cross(forceTerm[dB]);
deltaR[dBxAD] = forceTerm[dB].cross(deltaR[AD]);
if (FABS(pTrT2) > 1.0e-08) {
RealOpenMM delta2 = delta*two;
RealOpenMM delta2 = delta*2;
AmoebaReferenceForce::getCrossProduct(deltaR[BD], deltaR[CD], deltaR[BDxCD]);
AmoebaReferenceForce::getCrossProduct(deltaR[T], deltaR[CD], deltaR[TxCD]);
AmoebaReferenceForce::getCrossProduct(deltaR[AD], deltaR[BD], deltaR[ADxBD]);
AmoebaReferenceForce::getCrossProduct(deltaR[AD], deltaR[T], deltaR[ADxT]);
for (int ii = 0; ii < 3; ii++) {
RealOpenMM term = deltaR[BDxCD][ii] + delta2*deltaR[TxCD][ii];
forceTerm[dA][ii] += delta*deltaR[CDxdB][ii] + term*pTrT2;
term = deltaR[ADxBD][ii] + delta2*deltaR[ADxT][ii];
forceTerm[dC][ii] += delta*deltaR[dBxAD][ii] + term*pTrT2;
forceTerm[dD][ii] = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
}
deltaR[BDxCD] = forceTerm[dB].cross(deltaR[CD]);
deltaR[TxCD] = forceTerm[T].cross(deltaR[CD]);
deltaR[ADxBD] = forceTerm[AD].cross(deltaR[BD]);
deltaR[ADxT] = forceTerm[AD].cross(deltaR[T]);
RealVec term = deltaR[BDxCD] + deltaR[TxCD]*delta2;
forceTerm[dA] += deltaR[CDxdB]*delta + term*pTrT2;
term = deltaR[ADxBD] + deltaR[ADxT]*delta2;
forceTerm[dC] += deltaR[dBxAD]*delta + term*pTrT2;
forceTerm[dD] = -(forceTerm[dA] + forceTerm[dB] + forceTerm[dC]);
} else {
for (int ii = 0; ii < 3; ii++) {
forceTerm[dA][ii] += delta*deltaR[CDxdB][ii];
forceTerm[dC][ii] += delta*deltaR[dBxAD][ii];
forceTerm[dD][ii] = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
}
forceTerm[dA] += deltaR[CDxdB]*delta;
forceTerm[dC] += deltaR[dBxAD]*delta;
forceTerm[dD] = -(forceTerm[dA] + forceTerm[dB] + forceTerm[dC]);
}
// accumulate forces
for (int jj = 0; jj < 4; jj++) {
forces[jj][0] = forceTerm[jj][0];
forces[jj][1] = forceTerm[jj][1];
forces[jj][2] = forceTerm[jj][2];
}
for (int jj = 0; jj < 4; jj++)
forces[jj] = forceTerm[jj];
return energy;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceInPlaneAngleForce.h
View file @ 94fbbe9c
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -40,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceInPlaneAngleForce() {};
AmoebaReferenceInPlaneAngleForce() : usePeriodic(false) {};
/**---------------------------------------------------------------------------------------
......@@ -50,6 +50,16 @@ public:
~AmoebaReferenceInPlaneAngleForce() {};
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions.
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
Calculate Amoeba in-plane angle ixns (force and energy)
......@@ -88,6 +98,9 @@ public:
private:
bool usePeriodic;
RealVec boxVectors[3];
/**---------------------------------------------------------------------------------------
Get dEdT and energy prefactor given cosine of angle :: the calculation for different
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceOutOfPlaneBendForce.cpp
View file @ 94fbbe9c
......@@ -28,6 +28,13 @@
using std::vector;
using namespace OpenMM;
void AmoebaReferenceOutOfPlaneBendForce::setPeriodic(OpenMM::RealVec* vectors) {
usePeriodic = true;
boxVectors[0] = vectors[0];
boxVectors[1] = vectors[1];
boxVectors[2] = vectors[2];
}
/**---------------------------------------------------------------------------------------
Calculate Amoeba Out-Of-Plane-Bend ixn (force and energy)
......@@ -78,12 +85,21 @@ RealOpenMM AmoebaReferenceOutOfPlaneBendForce::calculateOutOfPlaneBendIxn(const
std::vector<RealOpenMM> deltaR[LastDeltaIndex];
for (int ii = 0; ii < LastDeltaIndex; ii++) {
deltaR[ii].resize(3);
}
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomA, deltaR[AB]);
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomC, deltaR[CB]);
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomD, deltaR[DB]);
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomA, deltaR[AD]);
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomC, deltaR[CD]);
}
if (usePeriodic) {
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomB, positionAtomA, deltaR[AB], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomB, positionAtomC, deltaR[CB], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomB, positionAtomD, deltaR[DB], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomD, positionAtomA, deltaR[AD], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomD, positionAtomC, deltaR[CD], boxVectors);
}
else {
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomA, deltaR[AB]);
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomC, deltaR[CB]);
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomD, deltaR[DB]);
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomA, deltaR[AD]);
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomC, deltaR[CD]);
}
RealOpenMM rDB2 = AmoebaReferenceForce::getNormSquared3(deltaR[DB]);
RealOpenMM rAD2 = AmoebaReferenceForce::getNormSquared3(deltaR[AD]);
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceOutOfPlaneBendForce.h
View file @ 94fbbe9c
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -40,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceOutOfPlaneBendForce() {};
AmoebaReferenceOutOfPlaneBendForce() : usePeriodic(false) {};
/**---------------------------------------------------------------------------------------
......@@ -50,6 +50,16 @@ public:
~AmoebaReferenceOutOfPlaneBendForce() {};
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions.
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
Calculate Amoeba out-of-plane-bend angle (force and energy)
......@@ -86,6 +96,9 @@ public:
private:
bool usePeriodic;
RealVec boxVectors[3];
/**---------------------------------------------------------------------------------------
Calculate Amoeba Out-Of-Plane-Bend ixn (force and energy)
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferencePiTorsionForce.cpp
View file @ 94fbbe9c
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -29,6 +28,13 @@
using std::vector;
using namespace OpenMM;
void AmoebaReferencePiTorsionForce::setPeriodic(OpenMM::RealVec* vectors) {
usePeriodic = true;
boxVectors[0] = vectors[0];
boxVectors[1] = vectors[1];
boxVectors[2] = vectors[2];
}
/**---------------------------------------------------------------------------------------
Calculate Amoeba pi-torsion ixn (force and energy)
......@@ -66,11 +72,19 @@ RealOpenMM AmoebaReferencePiTorsionForce::calculatePiTorsionIxn(const RealVec& p
std::vector<RealOpenMM> deltaR[LastDeltaIndex];
for (unsigned int ii = 0; ii < LastDeltaIndex; ii++) {
deltaR[ii].resize(3);
}
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomA, deltaR[AD]);
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomB, deltaR[BD]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomE, deltaR[EC]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomF, deltaR[FC]);
}
if (usePeriodic) {
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomD, positionAtomA, deltaR[AD], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomD, positionAtomB, deltaR[BD], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomC, positionAtomE, deltaR[EC], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomC, positionAtomF, deltaR[FC], boxVectors);
}
else {
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomA, deltaR[AD]);
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomB, deltaR[BD]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomE, deltaR[EC]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomF, deltaR[FC]);
}
enum { A, B, C, D, E, F, LastAtomIndex };
std::vector<RealOpenMM> d[LastAtomIndex];
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferencePiTorsionForce.h
View file @ 94fbbe9c
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -40,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferencePiTorsionForce() {};
AmoebaReferencePiTorsionForce() : usePeriodic(false) {};
/**---------------------------------------------------------------------------------------
......@@ -50,6 +50,16 @@ public:
~AmoebaReferencePiTorsionForce() {};
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions.
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
Calculate Amoeba torsion ixns (force and energy)
......@@ -85,6 +95,9 @@ public:
private:
bool usePeriodic;
RealVec boxVectors[3];
/**---------------------------------------------------------------------------------------
Calculate Amoeba pi-torsion ixn (force and energy)
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceStretchBendForce.cpp
View file @ 94fbbe9c
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -29,6 +29,13 @@
using std::vector;
using namespace OpenMM;
void AmoebaReferenceStretchBendForce::setPeriodic(OpenMM::RealVec* vectors) {
usePeriodic = true;
boxVectors[0] = vectors[0];
boxVectors[1] = vectors[1];
boxVectors[2] = vectors[2];
}
/**---------------------------------------------------------------------------------------
Calculate Amoeba stretch bend angle ixn (force and energy)
......@@ -76,8 +83,14 @@ RealOpenMM AmoebaReferenceStretchBendForce::calculateStretchBendIxn(const RealVe
for (unsigned int ii = 0; ii < LastDeltaIndex; ii++) {
deltaR[ii].resize(3);
}
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomA, deltaR[AB]);
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomC, deltaR[CB]);
if (usePeriodic) {
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomB, positionAtomA, deltaR[AB], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomB, positionAtomC, deltaR[CB], boxVectors);
}
else {
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomA, deltaR[AB]);
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomC, deltaR[CB]);
}
RealOpenMM rAB2 = AmoebaReferenceForce::getNormSquared3(deltaR[AB]);
RealOpenMM rAB = SQRT(rAB2);
RealOpenMM rCB2 = AmoebaReferenceForce::getNormSquared3(deltaR[CB]);
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceStretchBendForce.h
View file @ 94fbbe9c
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -40,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceStretchBendForce() {};
AmoebaReferenceStretchBendForce() : usePeriodic(false) {};
/**---------------------------------------------------------------------------------------
......@@ -50,6 +50,16 @@ public:
~AmoebaReferenceStretchBendForce() {};
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions.
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
Calculate Amoeba stretch bend ixns (force and energy)
......@@ -84,6 +94,9 @@ public:
private:
bool usePeriodic;
RealVec boxVectors[3];
/**---------------------------------------------------------------------------------------
Calculate Amoeba stretch bend angle ixn (force and energy)
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionTorsionForce.cpp
View file @ 94fbbe9c
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -28,6 +27,13 @@
using std::vector;
using namespace OpenMM;
void AmoebaReferenceTorsionTorsionForce::setPeriodic(OpenMM::RealVec* vectors) {
usePeriodic = true;
boxVectors[0] = vectors[0];
boxVectors[1] = vectors[1];
boxVectors[2] = vectors[2];
}
/**---------------------------------------------------------------------------------------
Load grid values from rectenclosing angles
......@@ -294,9 +300,16 @@ int AmoebaReferenceTorsionTorsionForce::checkTorsionSign(
deltaR[ii].resize(3);
}
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomA, deltaR[CA]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomB, deltaR[CB]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomD, deltaR[CD]);
if (usePeriodic) {
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomC, positionAtomA, deltaR[CA], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomC, positionAtomB, deltaR[CB], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomC, positionAtomD, deltaR[CD], boxVectors);
}
else {
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomA, deltaR[CA]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomB, deltaR[CB]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomD, deltaR[CD]);
}
RealOpenMM volume = deltaR[CA][0]*(deltaR[CB][1]*deltaR[CD][2] - deltaR[CB][2]*deltaR[CD][1]) +
deltaR[CB][0]*(deltaR[CD][1]*deltaR[CA][2] - deltaR[CD][2]*deltaR[CA][1]) +
......@@ -351,13 +364,24 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn(const
deltaR[ii].resize(3);
}
AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomB, deltaR[BA]);
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomC, deltaR[CB]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomD, deltaR[DC]);
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomE, deltaR[ED]);
AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomC, deltaR[CA]);
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomD, deltaR[DB]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomE, deltaR[EC]);
if (usePeriodic) {
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomA, positionAtomB, deltaR[BA], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomB, positionAtomC, deltaR[CB], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomC, positionAtomD, deltaR[DC], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomD, positionAtomE, deltaR[ED], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomA, positionAtomC, deltaR[CA], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomB, positionAtomD, deltaR[DB], boxVectors);
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomC, positionAtomE, deltaR[EC], boxVectors);
}
else {
AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomB, deltaR[BA]);
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomC, deltaR[CB]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomD, deltaR[DC]);
AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomE, deltaR[ED]);
AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomC, deltaR[CA]);
AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomD, deltaR[DB]);
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomE, deltaR[EC]);
}
std::vector<RealOpenMM> d[LastAtomIndex];
for (unsigned int ii = 0; ii < LastAtomIndex; ii++) {
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionTorsionForce.h
View file @ 94fbbe9c
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -40,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceTorsionTorsionForce() {};
AmoebaReferenceTorsionTorsionForce() : usePeriodic(false) {};
/**---------------------------------------------------------------------------------------
......@@ -50,6 +50,16 @@ public:
~AmoebaReferenceTorsionTorsionForce() {};
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions.
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic(OpenMM::RealVec* vectors);
/**---------------------------------------------------------------------------------------
Calculate Amoeba torsion-torsion ixns (force and energy)
......@@ -84,6 +94,9 @@ public:
private:
bool usePeriodic;
RealVec boxVectors[3];
/**---------------------------------------------------------------------------------------
Load grid values from rectangle enclosing angles
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaAngleForce.cpp
View file @ 94fbbe9c
......@@ -35,6 +35,7 @@
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/CustomAngleForce.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
......@@ -274,12 +275,63 @@ void testOneAngle() {
compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions, then compare to an identical custom force.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numParticles = 3;
for (int ii = 0; ii < numParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaAngleForce* amoebaAngleForce = new AmoebaAngleForce();
double angle = 100.0;
double quadraticK = 1.0;
double cubicK = 1.0e-01;
double quarticK = 1.0e-02;
double penticK = 1.0e-03;
double sexticK = 1.0e-04;
amoebaAngleForce->addAngle(0, 1, 2, angle, quadraticK);
amoebaAngleForce->setAmoebaGlobalAngleCubic(cubicK);
amoebaAngleForce->setAmoebaGlobalAngleQuartic(quarticK);
amoebaAngleForce->setAmoebaGlobalAnglePentic(penticK);
amoebaAngleForce->setAmoebaGlobalAngleSextic(sexticK);
amoebaAngleForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaAngleForce);
CustomAngleForce* customForce = new CustomAngleForce("k2*delta^2 + k3*delta^3 + k4*delta^4 + k5*delta^5 + k6*delta^6; delta=theta-theta0");
customForce->addGlobalParameter("theta0", angle*M_PI/180);
customForce->addGlobalParameter("k2", quadraticK*pow(180/M_PI, 2.0));
customForce->addGlobalParameter("k3", cubicK*pow(180/M_PI, 3.0));
customForce->addGlobalParameter("k4", quarticK*pow(180/M_PI, 4.0));
customForce->addGlobalParameter("k5", penticK*pow(180/M_PI, 5.0));
customForce->addGlobalParameter("k6", sexticK*pow(180/M_PI, 6.0));
customForce->addAngle(0, 1, 2);
customForce->setUsesPeriodicBoundaryConditions(true);
customForce->setForceGroup(1);
system.addForce(customForce);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0, 0, 2);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy, true, 1);
State s2 = context.getState(State::Forces | State::Energy, true, 2);
ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], 1e-5);
}
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestCudaAmoebaAngleForce running test..." << std::endl;
std::cout << "TestReferenceAmoebaAngleForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories();
testOneAngle();
testPeriodic();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaBondForce.cpp
View file @ 94fbbe9c
......@@ -35,6 +35,7 @@
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "openmm/CustomBondForce.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
......@@ -201,6 +202,49 @@ void testTwoBond() {
compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond");
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions, then compare to an identical custom force.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numParticles = 2;
for (int ii = 0; ii < numParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaBondForce* amoebaBondForce = new AmoebaBondForce();
double bondLength = 1.5;
double quadraticK = 1.0;
double cubicK = 2.0;
double quarticK = 3.0;
amoebaBondForce->setAmoebaGlobalBondCubic(cubicK);
amoebaBondForce->setAmoebaGlobalBondQuartic(quarticK);
amoebaBondForce->addBond(0, 1, bondLength, quadraticK);
amoebaBondForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaBondForce);
CustomBondForce* customForce = new CustomBondForce("k2*delta^2 + k3*delta^3 + k4*delta^4; delta=r-r0");
customForce->addGlobalParameter("r0", bondLength);
customForce->addGlobalParameter("k2", quadraticK);
customForce->addGlobalParameter("k3", cubicK);
customForce->addGlobalParameter("k4", quarticK);
customForce->addBond(0, 1);
customForce->setUsesPeriodicBoundaryConditions(true);
customForce->setForceGroup(1);
system.addForce(customForce);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numParticles);
positions[0] = Vec3(0, 2, 0);
positions[1] = Vec3(0, 0, 0);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy, true, 1);
State s2 = context.getState(State::Forces | State::Energy, true, 2);
ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], 1e-5);
}
int main(int numberOfArguments, char* argv[]) {
try {
......@@ -208,6 +252,7 @@ int main(int numberOfArguments, char* argv[]) {
registerAmoebaReferenceKernelFactories();
//testOneBond();
testTwoBond();
testPeriodic();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaInPlaneAngleForce.cpp
View file @ 94fbbe9c
......@@ -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 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -354,12 +354,57 @@ void testOneAngle() {
compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numberOfParticles = 4;
for (int ii = 0; ii < numberOfParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaInPlaneAngleForce* amoebaInPlaneAngleForce = new AmoebaInPlaneAngleForce();
double angle = 65.0;
double quadraticK = 1.0;
double cubicK = 0.0e-01;
double quarticK = 0.0e-02;
double penticK = 0.0e-03;
double sexticK = 0.0e-04;
amoebaInPlaneAngleForce->addAngle(0, 1, 2, 3, angle, quadraticK);
amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleCubic(cubicK);
amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleQuartic(quarticK);
amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAnglePentic(penticK);
amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleSextic(sexticK);
amoebaInPlaneAngleForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaInPlaneAngleForce);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0, 0, 1);
positions[3] = Vec3(1, 1, 1);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy);
// Move one atom to a position that should give identical results.
positions[2] = Vec3(0, 0, -2);
context.setPositions(positions);
State s2 = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestReferenceAmoebaInPlaneAngleForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories();
testOneAngle();
testPeriodic();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaOutOfPlaneBendForce.cpp
View file @ 94fbbe9c
......@@ -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 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -461,12 +461,50 @@ void testOneOutOfPlaneBend2(int setId) {
}
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numberOfParticles = 4;
for (int ii = 0; ii < numberOfParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce();
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic(0.5600000E-04);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic(-0.7000000E-06);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07);
double kOutOfPlaneBend = 0.328682196E-01;
amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend);
amoebaOutOfPlaneBendForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaOutOfPlaneBendForce);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(1, 0, 0);
positions[2] = Vec3(0, 1, 0);
positions[3] = Vec3(0, 0, 1);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy);
// Move one atom to a position that should give identical results.
positions[3] = Vec3(0, 0, -2);
context.setPositions(positions);
State s2 = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestReferenceAmoebaOutOfPlaneBendForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories();
testOneOutOfPlaneBend();
testPeriodic();
//testOneOutOfPlaneBend2(atoi(argv[1]));
//for (int ii = 1; ii <= 6; ii++) {
// testOneOutOfPlaneBend2(ii);
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaPiTorsionForce.cpp
View file @ 94fbbe9c
......@@ -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 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -290,12 +290,48 @@ void testOnePiTorsion() {
compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion");
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numberOfParticles = 6;
for (int ii = 0; ii < numberOfParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaPiTorsionForce* amoebaPiTorsionForce = new AmoebaPiTorsionForce();
double kTorsion = 6.85;
amoebaPiTorsionForce->addPiTorsion(0, 1, 2, 3, 4, 5, kTorsion);
amoebaPiTorsionForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaPiTorsionForce);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0, 0, 0.5);
positions[3] = Vec3(0.4, 0.4, 0.4);
positions[4] = Vec3(1, 0, 1);
positions[5] = Vec3(1, 1, 0);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy);
// Move one atom to a position that should give identical results.
positions[0] = Vec3(0, -2, 0);
context.setPositions(positions);
State s2 = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestReferenceAmoebaPiTorsionForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories();
testOnePiTorsion();
testPeriodic();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaStretchBendForce.cpp
View file @ 94fbbe9c
......@@ -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 Stanford University and the Authors. *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -270,12 +270,48 @@ void testOneStretchBend() {
compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
}
void testPeriodic() {
// Create a force that uses periodic boundary conditions.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
int numberOfParticles = 3;
for (int ii = 0; ii < numberOfParticles; ii++)
system.addParticle(1.0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
AmoebaStretchBendForce* amoebaStretchBendForce = new AmoebaStretchBendForce();
double abLength = 0.144800000E+01;
double cbLength = 0.101500000E+01;
double angleStretchBend = 0.108500000E+03*DegreesToRadians;
double kStretchBend = 1.0;
amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend);
amoebaStretchBendForce->setUsesPeriodicBoundaryConditions(true);
system.addForce(amoebaStretchBendForce);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
positions[2] = Vec3(0, 0, 1);
context.setPositions(positions);
State s1 = context.getState(State::Forces | State::Energy);
// Move one atom to a position that should give identical results.
positions[2] = Vec3(0, 0, -2);
context.setPositions(positions);
State s2 = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
for (int i = 0; i < numberOfParticles; i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestReferenceAmoebaStretchBendForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories();
testOneStretchBend();
testPeriodic();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
......
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