Skip to content

GitLab

Commit 41cd79a5 authored by peastman's avatar peastman
Browse files

Cleaned up formatting of AMOEBA reference code

parent 25a308e6
Hide whitespace changes
Inline Side-by-side
Showing with 4541 additions and 4541 deletions
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionTorsionForce.cpp
View file @ 41cd79a5
......@@ -50,7 +50,7 @@ using namespace OpenMM;
void AmoebaReferenceTorsionTorsionForce::loadGridValuesFromEnclosingRectangle(
const std::vector< std::vector< std::vector<RealOpenMM> > >& grid,
RealOpenMM angle1, RealOpenMM angle2, RealOpenMM corners[2][2],
RealOpenMM* fValues, RealOpenMM* fValues1, RealOpenMM* fValues2, RealOpenMM* fValues12 ) const {
RealOpenMM* fValues, RealOpenMM* fValues1, RealOpenMM* fValues2, RealOpenMM* fValues12) const {
// ---------------------------------------------------------------------------------------
......@@ -62,8 +62,8 @@ void AmoebaReferenceTorsionTorsionForce::loadGridValuesFromEnclosingRectangle(
unsigned int gridSize = grid.size();
double gridSpacingI = static_cast<double>(gridSize-1)/360.0;
const int X_gridIndex = (int) ( (angle1 - grid[0][0][0] )*gridSpacingI + 1.0e-06);
const int Y_gridIndex = (int) ( (angle2 - grid[0][0][1] )*gridSpacingI + 1.0e-06);
const int X_gridIndex = (int) ((angle1 - grid[0][0][0])*gridSpacingI + 1.0e-06);
const int Y_gridIndex = (int) ((angle2 - grid[0][0][1])*gridSpacingI + 1.0e-06);
// get coordinates of corner indices
......@@ -73,15 +73,15 @@ void AmoebaReferenceTorsionTorsionForce::loadGridValuesFromEnclosingRectangle(
corners[1][0] = grid[X_gridIndex][Y_gridIndex ][1];
corners[1][1] = grid[X_gridIndex+1][Y_gridIndex+1][1];
for( int ii = 0; ii < 4; ii++ ){
for (int ii = 0; ii < 4; ii++) {
int gridX = X_gridIndex;
int gridY = Y_gridIndex;
if( ii == 1 ){
if (ii == 1) {
gridX++;
} else if( ii == 2 ){
} else if (ii == 2) {
gridX++;
gridY++;
} else if( ii == 3 ){
} else if (ii == 3) {
gridY++;
}
......@@ -117,9 +117,9 @@ void AmoebaReferenceTorsionTorsionForce::loadGridValuesFromEnclosingRectangle(
--------------------------------------------------------------------------------------- */
void AmoebaReferenceTorsionTorsionForce::getBicubicCoefficientMatrix( const RealOpenMM* y,
void AmoebaReferenceTorsionTorsionForce::getBicubicCoefficientMatrix(const RealOpenMM* y,
const RealOpenMM* y1, const RealOpenMM* y2, const RealOpenMM* y12, const RealOpenMM d1, const RealOpenMM d2,
RealOpenMM c[4][4] ) const {
RealOpenMM c[4][4]) const {
// ---------------------------------------------------------------------------------------
......@@ -159,9 +159,9 @@ void AmoebaReferenceTorsionTorsionForce::getBicubicCoefficientMatrix( const Real
// pack y, y1, y2, y12 into single vector of dimension 16
std::vector<RealOpenMM> x( 16 );
std::vector<RealOpenMM> x(16);
RealOpenMM d1d2 = d1*d2;
for( int ii = 0; ii < 4; ii++ ){
for (int ii = 0; ii < 4; ii++) {
x[ii] = y[ii];
x[ii+4] = y1[ii]*d1;
x[ii+8] = y2[ii]*d2;
......@@ -172,13 +172,13 @@ void AmoebaReferenceTorsionTorsionForce::getBicubicCoefficientMatrix( const Real
int rowIndex = 0;
int colIndex = 0;
for( int ii = 0; ii < 16; ii++ ){
for (int ii = 0; ii < 16; ii++) {
RealOpenMM sum = weightMatrix[0][ii]*x[0];
for( int jj = 1; jj < 16; jj++ ){
for (int jj = 1; jj < 16; jj++) {
sum += weightMatrix[jj][ii]*x[jj];
}
c[rowIndex][colIndex++] = sum;
if( !(colIndex % 4) ){
if (!(colIndex % 4)) {
colIndex = 0;
rowIndex++;
}
......@@ -216,12 +216,12 @@ void AmoebaReferenceTorsionTorsionForce::getBicubicCoefficientMatrix( const Real
--------------------------------------------------------------------------------------- */
void AmoebaReferenceTorsionTorsionForce::getBicubicValues(
void AmoebaReferenceTorsionTorsionForce::getBicubicValues(
const RealOpenMM* y, const RealOpenMM* y1, const RealOpenMM* y2, const RealOpenMM* y12,
const RealOpenMM x1Lower, const RealOpenMM x1Upper,
const RealOpenMM x2Lower, const RealOpenMM x2Upper,
const RealOpenMM gridValue1, const RealOpenMM gridValue2,
RealOpenMM* functionValue, RealOpenMM* functionValue1, RealOpenMM* functionValue2 ) const {
RealOpenMM* functionValue, RealOpenMM* functionValue1, RealOpenMM* functionValue2) const {
// ---------------------------------------------------------------------------------------
......@@ -236,7 +236,7 @@ void AmoebaReferenceTorsionTorsionForce::getBicubicValues(
// get coefficent matrix
RealOpenMM coefficientMatrix[4][4];
getBicubicCoefficientMatrix( y, y1, y2, y12, x1Upper-x1Lower, x2Upper-x2Lower, coefficientMatrix );
getBicubicCoefficientMatrix(y, y1, y2, y12, x1Upper-x1Lower, x2Upper-x2Lower, coefficientMatrix);
// apply coefficent matrix
......@@ -247,10 +247,10 @@ void AmoebaReferenceTorsionTorsionForce::getBicubicValues(
*functionValue1 = zero;
*functionValue2 = zero;
for( int ii = 3; ii >= 0; ii-- ){
*functionValue = t*(*functionValue) + ( ( coefficientMatrix[ii][3]*u + coefficientMatrix[ii][2] )*u + coefficientMatrix[ii][1])*u + coefficientMatrix[ii][0];
*functionValue1 = u*(*functionValue1) + ( three*coefficientMatrix[3][ii]*t + two*coefficientMatrix[2][ii] )*t + coefficientMatrix[1][ii];
*functionValue2 = t*(*functionValue2) + ( three*coefficientMatrix[ii][3]*u + two*coefficientMatrix[ii][2] )*u + coefficientMatrix[ii][1];
for (int ii = 3; ii >= 0; ii--) {
*functionValue = t*(*functionValue) + ( (coefficientMatrix[ii][3]*u + coefficientMatrix[ii][2])*u + coefficientMatrix[ii][1])*u + coefficientMatrix[ii][0];
*functionValue1 = u*(*functionValue1) + (three*coefficientMatrix[3][ii]*t + two*coefficientMatrix[2][ii])*t + coefficientMatrix[1][ii];
*functionValue2 = t*(*functionValue2) + (three*coefficientMatrix[ii][3]*u + two*coefficientMatrix[ii][2])*u + coefficientMatrix[ii][1];
}
*functionValue1 /= (x1Upper - x1Lower);
......@@ -273,9 +273,9 @@ void AmoebaReferenceTorsionTorsionForce::getBicubicValues(
--------------------------------------------------------------------------------------- */
int AmoebaReferenceTorsionTorsionForce::checkTorsionSign(
int AmoebaReferenceTorsionTorsionForce::checkTorsionSign(
const RealVec& positionAtomA, const RealVec& positionAtomB,
const RealVec& positionAtomC, const RealVec& positionAtomD ) const {
const RealVec& positionAtomC, const RealVec& positionAtomD) const {
// ---------------------------------------------------------------------------------------
......@@ -290,17 +290,17 @@ int AmoebaReferenceTorsionTorsionForce::checkTorsionSign(
enum { CA, CB, CD, LastDeltaIndex };
std::vector<RealOpenMM> deltaR[LastDeltaIndex];
for( unsigned int ii = 0; ii < LastDeltaIndex; ii++ ){
for (unsigned int ii = 0; ii < LastDeltaIndex; ii++) {
deltaR[ii].resize(3);
}
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomA, deltaR[CA] );
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomB, deltaR[CB] );
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomD, deltaR[CD] );
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] ) +
deltaR[CD][0]*( deltaR[CA][1]*deltaR[CB][2] - deltaR[CA][2]*deltaR[CB][1] );
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]) +
deltaR[CD][0]*(deltaR[CA][1]*deltaR[CB][2] - deltaR[CA][2]*deltaR[CB][1]);
return (volume >= zero ? one : -one);
......@@ -324,11 +324,11 @@ int AmoebaReferenceTorsionTorsionForce::checkTorsionSign(
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( const RealVec& positionAtomA, const RealVec& positionAtomB,
const RealVec& positionAtomC, const RealVec& positionAtomD,
const RealVec& positionAtomE, const RealVec* positionChiralCheckAtom,
const std::vector< std::vector< std::vector<RealOpenMM> > >& grid,
RealVec* forces ) const {
RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn(const RealVec& positionAtomA, const RealVec& positionAtomB,
const RealVec& positionAtomC, const RealVec& positionAtomD,
const RealVec& positionAtomE, const RealVec* positionChiralCheckAtom,
const std::vector< std::vector< std::vector<RealOpenMM> > >& grid,
RealVec* forces) const {
// ---------------------------------------------------------------------------------------
......@@ -347,83 +347,83 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( const
enum { BA, CB, DC, ED, T, U, V, UxV, CA, DB, EC, dT, dU, dU2, dV2, LastDeltaIndex };
std::vector<RealOpenMM> deltaR[LastDeltaIndex];
for( unsigned int ii = 0; ii < LastDeltaIndex; ii++ ){
for (unsigned int ii = 0; ii < LastDeltaIndex; ii++) {
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] );
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++ ){
for (unsigned int ii = 0; ii < LastAtomIndex; ii++) {
d[ii].resize(3);
}
AmoebaReferenceForce::getCrossProduct( deltaR[BA], deltaR[CB], deltaR[T] );
AmoebaReferenceForce::getCrossProduct( deltaR[CB], deltaR[DC], deltaR[U] );
AmoebaReferenceForce::getCrossProduct( deltaR[DC], deltaR[ED], deltaR[V] );
AmoebaReferenceForce::getCrossProduct(deltaR[BA], deltaR[CB], deltaR[T]);
AmoebaReferenceForce::getCrossProduct(deltaR[CB], deltaR[DC], deltaR[U]);
AmoebaReferenceForce::getCrossProduct(deltaR[DC], deltaR[ED], deltaR[V]);
AmoebaReferenceForce::getCrossProduct( deltaR[U], deltaR[V], deltaR[UxV] );
AmoebaReferenceForce::getCrossProduct(deltaR[U], deltaR[V], deltaR[UxV]);
RealOpenMM rT2 = AmoebaReferenceForce::getNormSquared3( deltaR[T] );
RealOpenMM rU2 = AmoebaReferenceForce::getNormSquared3( deltaR[U] );
RealOpenMM rV2 = AmoebaReferenceForce::getNormSquared3( deltaR[V] );
RealOpenMM rUrV = SQRT( rU2*rV2 );
RealOpenMM rT2 = AmoebaReferenceForce::getNormSquared3(deltaR[T]);
RealOpenMM rU2 = AmoebaReferenceForce::getNormSquared3(deltaR[U]);
RealOpenMM rV2 = AmoebaReferenceForce::getNormSquared3(deltaR[V]);
RealOpenMM rUrV = SQRT(rU2*rV2);
RealOpenMM rTrU = SQRT( rT2*rU2 );
RealOpenMM rTrU = SQRT(rT2*rU2);
if( rTrU <= zero || rUrV <= zero ){
if (rTrU <= zero || rUrV <= zero) {
return zero;
}
RealOpenMM rCB = AmoebaReferenceForce::getNorm3( deltaR[CB] );
RealOpenMM cosine1 = AmoebaReferenceForce::getDotProduct3( deltaR[T], deltaR[U] );
RealOpenMM rCB = AmoebaReferenceForce::getNorm3(deltaR[CB]);
RealOpenMM cosine1 = AmoebaReferenceForce::getDotProduct3(deltaR[T], deltaR[U]);
cosine1 /= rTrU;
RealOpenMM angle1;
if( cosine1 <= -one ){
if (cosine1 <= -one) {
angle1 = PI_M*RADIAN;
} else if( cosine1 >= one ){
} else if (cosine1 >= one) {
angle1 = zero;
} else {
angle1 = RADIAN*ACOS(cosine1);
}
RealOpenMM sign = AmoebaReferenceForce::getDotProduct3( deltaR[BA], deltaR[U] );
if( sign < zero ){
RealOpenMM sign = AmoebaReferenceForce::getDotProduct3(deltaR[BA], deltaR[U]);
if (sign < zero) {
angle1 = -angle1;
}
// value1 = angle1;
RealOpenMM rDC = AmoebaReferenceForce::getNorm3( deltaR[DC] );
RealOpenMM cosine2 = AmoebaReferenceForce::getDotProduct3( deltaR[U], deltaR[V] );
RealOpenMM rDC = AmoebaReferenceForce::getNorm3(deltaR[DC]);
RealOpenMM cosine2 = AmoebaReferenceForce::getDotProduct3(deltaR[U], deltaR[V]);
cosine2 /= rUrV;
RealOpenMM angle2;
if( cosine2 <= -one ){
if (cosine2 <= -one) {
angle2 = PI_M*RADIAN;
} else if( cosine1 >= one ){
} else if (cosine1 >= one) {
angle2 = zero;
} else {
angle2 = RADIAN*ACOS(cosine2);
}
sign = AmoebaReferenceForce::getDotProduct3( deltaR[CB], deltaR[V] );
if( sign < zero ){
sign = AmoebaReferenceForce::getDotProduct3(deltaR[CB], deltaR[V]);
if (sign < zero) {
angle2 = -angle2;
}
// swap signs of angles if chirality at central atom
// is 'negative'
if( positionChiralCheckAtom ){
sign = checkTorsionSign( *positionChiralCheckAtom, positionAtomB, positionAtomC, positionAtomD );
if( sign < zero ){
if (positionChiralCheckAtom) {
sign = checkTorsionSign(*positionChiralCheckAtom, positionAtomB, positionAtomC, positionAtomD);
if (sign < zero) {
angle1 = -angle1;
angle2 = -angle2;
}
......@@ -436,7 +436,7 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( const
RealOpenMM corners[2][2];
RealOpenMM eValues[4][4];
enum { E0, E1, E2, E12, LastEIndex };
loadGridValuesFromEnclosingRectangle( grid, angle1, angle2, corners, eValues[E0], eValues[E1], eValues[E2], eValues[E12] );
loadGridValuesFromEnclosingRectangle(grid, angle1, angle2, corners, eValues[E0], eValues[E1], eValues[E2], eValues[E12]);
// get coordinates of point in grid closest to angle1 & angle2
......@@ -446,16 +446,16 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( const
RealOpenMM gridEnergy;
RealOpenMM dEdAngle1;
RealOpenMM dEdAngle2;
AmoebaReferenceTorsionTorsionForce::getBicubicValues(
AmoebaReferenceTorsionTorsionForce::getBicubicValues(
eValues[E0], eValues[E1], eValues[E2], eValues[E12],
corners[0][0], corners[0][1], corners[1][0], corners[1][1],
angle1, angle2, &gridEnergy, &dEdAngle1, &dEdAngle2 );
angle1, angle2, &gridEnergy, &dEdAngle1, &dEdAngle2);
dEdAngle1 = sign*RADIAN*dEdAngle1;
dEdAngle2 = sign*RADIAN*dEdAngle2;
AmoebaReferenceForce::getCrossProduct( deltaR[T], deltaR[CB], deltaR[dT] );
AmoebaReferenceForce::getCrossProduct( deltaR[U], deltaR[CB], deltaR[dU] );
AmoebaReferenceForce::getCrossProduct(deltaR[T], deltaR[CB], deltaR[dT]);
AmoebaReferenceForce::getCrossProduct(deltaR[U], deltaR[CB], deltaR[dU]);
RealOpenMM factorT = dEdAngle1/(rCB*rT2);
RealOpenMM factorU = -dEdAngle1/(rCB*rU2);
......@@ -468,18 +468,18 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( const
deltaR[dU][1] *= factorU;
deltaR[dU][2] *= factorU;
AmoebaReferenceForce::getCrossProduct( deltaR[dT], deltaR[CB], d[A] );
AmoebaReferenceForce::getCrossProduct( deltaR[dU], deltaR[CB], d[D] );
AmoebaReferenceForce::getCrossProduct(deltaR[dT], deltaR[CB], d[A]);
AmoebaReferenceForce::getCrossProduct(deltaR[dU], deltaR[CB], d[D]);
std::vector<RealOpenMM> tmp[3];
for( unsigned int ii = 0; ii < 3; ii++ ){
for (unsigned int ii = 0; ii < 3; ii++) {
tmp[ii].resize(3);
}
AmoebaReferenceForce::getCrossProduct( deltaR[CA], deltaR[dT], d[B] );
AmoebaReferenceForce::getCrossProduct( deltaR[dU], deltaR[DC], tmp[0] );
AmoebaReferenceForce::getCrossProduct(deltaR[CA], deltaR[dT], d[B]);
AmoebaReferenceForce::getCrossProduct(deltaR[dU], deltaR[DC], tmp[0]);
AmoebaReferenceForce::getCrossProduct( deltaR[dT], deltaR[BA], d[C] );
AmoebaReferenceForce::getCrossProduct( deltaR[DB], deltaR[dU], tmp[1]);
AmoebaReferenceForce::getCrossProduct(deltaR[dT], deltaR[BA], d[C] );
AmoebaReferenceForce::getCrossProduct(deltaR[DB], deltaR[dU], tmp[1]);
d[B][0] += tmp[0][0];
d[B][1] += tmp[0][1];
......@@ -491,8 +491,8 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( const
// angle2
AmoebaReferenceForce::getCrossProduct( deltaR[U], deltaR[DC], deltaR[dU2] );
AmoebaReferenceForce::getCrossProduct( deltaR[V], deltaR[DC], deltaR[dV2] );
AmoebaReferenceForce::getCrossProduct(deltaR[U], deltaR[DC], deltaR[dU2]);
AmoebaReferenceForce::getCrossProduct(deltaR[V], deltaR[DC], deltaR[dV2]);
RealOpenMM factorU2 = dEdAngle2/(rDC*rU2);
RealOpenMM factorV2 = -dEdAngle2/(rDC*rV2);
......@@ -507,12 +507,12 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( const
// dB2
AmoebaReferenceForce::getCrossProduct( deltaR[dU2], deltaR[DC], tmp[0] );
AmoebaReferenceForce::getCrossProduct(deltaR[dU2], deltaR[DC], tmp[0] );
// dC2
AmoebaReferenceForce::getCrossProduct( deltaR[DB], deltaR[dU2], tmp[1] );
AmoebaReferenceForce::getCrossProduct( deltaR[dV2], deltaR[ED], tmp[2] );
AmoebaReferenceForce::getCrossProduct(deltaR[DB], deltaR[dU2], tmp[1] );
AmoebaReferenceForce::getCrossProduct(deltaR[dV2], deltaR[ED], tmp[2] );
d[B][0] += tmp[0][0];
d[B][1] += tmp[0][1];
......@@ -524,8 +524,8 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( const
// dD2
AmoebaReferenceForce::getCrossProduct( deltaR[dU2], deltaR[CB], tmp[0] );
AmoebaReferenceForce::getCrossProduct( deltaR[EC], deltaR[dV2], tmp[1] );
AmoebaReferenceForce::getCrossProduct(deltaR[dU2], deltaR[CB], tmp[0] );
AmoebaReferenceForce::getCrossProduct(deltaR[EC], deltaR[dV2], tmp[1] );
d[D][0] += tmp[0][0] + tmp[1][0];
d[D][1] += tmp[0][1] + tmp[1][1];
......@@ -533,13 +533,13 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( const
// dE
AmoebaReferenceForce::getCrossProduct( deltaR[dV2], deltaR[DC], d[E] );
AmoebaReferenceForce::getCrossProduct(deltaR[dV2], deltaR[DC], d[E] );
// ---------------------------------------------------------------------------------------
// add in forces
for( int jj = 0; jj < LastAtomIndex; jj++ ){
for (int jj = 0; jj < LastAtomIndex; jj++) {
forces[jj][0] = d[jj][0];
forces[jj][1] = d[jj][1];
forces[jj][2] = d[jj][2];
......@@ -550,16 +550,16 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( const
return gridEnergy;
}
RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateForceAndEnergy( int numTorsionTorsions, vector<RealVec>& posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<int>& particle5,
const std::vector<int>& chiralCheckAtom,
const std::vector<int>& gridIndices,
const std::vector< std::vector< std::vector< std::vector<RealOpenMM> > > >& torsionTorsionGrids,
vector<RealVec>& forceData ) const {
RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateForceAndEnergy(int numTorsionTorsions, vector<RealVec>& posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<int>& particle5,
const std::vector<int>& chiralCheckAtom,
const std::vector<int>& gridIndices,
const std::vector< std::vector< std::vector< std::vector<RealOpenMM> > > >& torsionTorsionGrids,
vector<RealVec>& forceData) const {
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numTorsionTorsions); ii++) {
......@@ -575,19 +575,19 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateForceAndEnergy( int numT
RealVec forces[5];
RealVec* chiralCheckAtom;
if( chiralCheckAtomIndex > -1 ){
if (chiralCheckAtomIndex > -1) {
chiralCheckAtom = &posData[chiralCheckAtomIndex];
} else {
chiralCheckAtom = NULL;
}
energy += calculateTorsionTorsionIxn( posData[particle1Index], posData[particle2Index],
posData[particle3Index], posData[particle4Index],
posData[particle5Index], chiralCheckAtom, torsionTorsionGrids[gridIndex],
forces );
energy += calculateTorsionTorsionIxn(posData[particle1Index], posData[particle2Index],
posData[particle3Index], posData[particle4Index],
posData[particle5Index], chiralCheckAtom, torsionTorsionGrids[gridIndex],
forces);
// accumulate forces
for( int jj = 0; jj < 3; jj++ ){
for (int jj = 0; jj < 3; jj++) {
forceData[particle1Index][jj] -= forces[0][jj];
forceData[particle2Index][jj] -= forces[1][jj];
forceData[particle3Index][jj] -= forces[2][jj];
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionTorsionForce.h
View file @ 41cd79a5
......@@ -40,7 +40,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceTorsionTorsionForce( ){};
AmoebaReferenceTorsionTorsionForce() {};
/**---------------------------------------------------------------------------------------
......@@ -48,7 +48,7 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferenceTorsionTorsionForce( ){};
~AmoebaReferenceTorsionTorsionForce() {};
/**---------------------------------------------------------------------------------------
......@@ -71,16 +71,16 @@ public:
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numTorsionTorsions, std::vector<OpenMM::RealVec>& posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<int>& particle5,
const std::vector<int>& chiralCheckAtom,
const std::vector<int>& gridIndices,
const std::vector< std::vector< std::vector< std::vector<RealOpenMM> > > >& torsionTorsionGrids,
std::vector<OpenMM::RealVec>& forceData ) const;
RealOpenMM calculateForceAndEnergy(int numTorsionTorsions, std::vector<OpenMM::RealVec>& posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector<int>& particle5,
const std::vector<int>& chiralCheckAtom,
const std::vector<int>& gridIndices,
const std::vector< std::vector< std::vector< std::vector<RealOpenMM> > > >& torsionTorsionGrids,
std::vector<OpenMM::RealVec>& forceData) const;
private:
......@@ -105,7 +105,7 @@ private:
void loadGridValuesFromEnclosingRectangle(
const std::vector< std::vector< std::vector<RealOpenMM> > >& grid,
RealOpenMM angle1, RealOpenMM angle2, RealOpenMM corners[2][2],
RealOpenMM* fValues, RealOpenMM* fValues1, RealOpenMM* fValues2, RealOpenMM* fValues12 ) const;
RealOpenMM* fValues, RealOpenMM* fValues1, RealOpenMM* fValues2, RealOpenMM* fValues12) const;
/**---------------------------------------------------------------------------------------
......@@ -128,8 +128,8 @@ private:
--------------------------------------------------------------------------------------- */
void getBicubicCoefficientMatrix( const RealOpenMM* y, const RealOpenMM* y1, const RealOpenMM* y2, const RealOpenMM* y12,
const RealOpenMM d1, const RealOpenMM d2, RealOpenMM c[4][4] ) const;
void getBicubicCoefficientMatrix(const RealOpenMM* y, const RealOpenMM* y1, const RealOpenMM* y2, const RealOpenMM* y12,
const RealOpenMM d1, const RealOpenMM d2, RealOpenMM c[4][4]) const;
/**---------------------------------------------------------------------------------------
......@@ -167,7 +167,7 @@ private:
const RealOpenMM x1Lower, const RealOpenMM x1Upper,
const RealOpenMM x2Lower, const RealOpenMM x2Upper,
const RealOpenMM gridValue1, const RealOpenMM gridValue2,
RealOpenMM* functionValue, RealOpenMM* functionValue1, RealOpenMM* functionValue2 ) const;
RealOpenMM* functionValue, RealOpenMM* functionValue1, RealOpenMM* functionValue2) const;
/**---------------------------------------------------------------------------------------
......@@ -183,8 +183,8 @@ private:
--------------------------------------------------------------------------------------- */
int checkTorsionSign( const OpenMM::RealVec& positionAtomA, const OpenMM::RealVec& positionAtomB,
const OpenMM::RealVec& positionAtomC, const OpenMM::RealVec& positionAtomD ) const;
int checkTorsionSign(const OpenMM::RealVec& positionAtomA, const OpenMM::RealVec& positionAtomB,
const OpenMM::RealVec& positionAtomC, const OpenMM::RealVec& positionAtomD) const;
/**---------------------------------------------------------------------------------------
......@@ -204,11 +204,11 @@ private:
--------------------------------------------------------------------------------------- */
RealOpenMM calculateTorsionTorsionIxn( const OpenMM::RealVec& positionAtomA, const OpenMM::RealVec& positionAtomB,
const OpenMM::RealVec& positionAtomC, const OpenMM::RealVec& positionAtomD,
const OpenMM::RealVec& positionAtomE, const OpenMM::RealVec* chiralCheckAtom,
const std::vector< std::vector< std::vector<RealOpenMM> > >& grid,
OpenMM::RealVec* forces ) const;
RealOpenMM calculateTorsionTorsionIxn(const OpenMM::RealVec& positionAtomA, const OpenMM::RealVec& positionAtomB,
const OpenMM::RealVec& positionAtomC, const OpenMM::RealVec& positionAtomD,
const OpenMM::RealVec& positionAtomE, const OpenMM::RealVec* chiralCheckAtom,
const std::vector< std::vector< std::vector<RealOpenMM> > >& grid,
OpenMM::RealVec* forces) const;
};
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp
View file @ 41cd79a5
......@@ -31,32 +31,32 @@
using std::vector;
using namespace OpenMM;
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( ) : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9) {
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce() : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9) {
setTaperCoefficients( _cutoff );
setSigmaCombiningRule( "ARITHMETIC" );
setEpsilonCombiningRule( "GEOMETRIC" );
setTaperCoefficients(_cutoff);
setSigmaCombiningRule("ARITHMETIC");
setEpsilonCombiningRule("GEOMETRIC");
}
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule ) : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9) {
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce(const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule) : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9) {
setTaperCoefficients( _cutoff );
setSigmaCombiningRule( sigmaCombiningRule );
setEpsilonCombiningRule( epsilonCombiningRule );
setTaperCoefficients(_cutoff);
setSigmaCombiningRule(sigmaCombiningRule);
setEpsilonCombiningRule(epsilonCombiningRule);
}
AmoebaReferenceVdwForce::NonbondedMethod AmoebaReferenceVdwForce::getNonbondedMethod( void ) const {
AmoebaReferenceVdwForce::NonbondedMethod AmoebaReferenceVdwForce::getNonbondedMethod() const {
return _nonbondedMethod;
}
void AmoebaReferenceVdwForce::setNonbondedMethod( AmoebaReferenceVdwForce::NonbondedMethod nonbondedMethod ){
void AmoebaReferenceVdwForce::setNonbondedMethod(AmoebaReferenceVdwForce::NonbondedMethod nonbondedMethod) {
_nonbondedMethod = nonbondedMethod;
}
void AmoebaReferenceVdwForce::setTaperCoefficients( double cutoff ){
void AmoebaReferenceVdwForce::setTaperCoefficients(double cutoff) {
_taperCutoff = cutoff*_taperCutoffFactor;
if( _taperCutoff != cutoff ){
if (_taperCutoff != cutoff) {
_taperCoefficients[C3] = 10.0/pow(_taperCutoff - cutoff, 3.0);
_taperCoefficients[C4] = 15.0/pow(_taperCutoff - cutoff, 4.0);
_taperCoefficients[C5] = 6.0/pow(_taperCutoff - cutoff, 5.0);
......@@ -67,12 +67,12 @@ void AmoebaReferenceVdwForce::setTaperCoefficients( double cutoff ){
}
}
void AmoebaReferenceVdwForce::setCutoff( double cutoff ){
void AmoebaReferenceVdwForce::setCutoff(double cutoff) {
_cutoff = cutoff;
setTaperCoefficients( _cutoff );
setTaperCoefficients(_cutoff);
}
double AmoebaReferenceVdwForce::getCutoff( void ) const {
double AmoebaReferenceVdwForce::getCutoff() const {
return _cutoff;
}
......@@ -82,35 +82,35 @@ void AmoebaReferenceVdwForce::setPeriodicBox(OpenMM::RealVec* vectors) {
_periodicBoxVectors[2] = vectors[2];
}
void AmoebaReferenceVdwForce::setSigmaCombiningRule( const std::string& sigmaCombiningRule ){
void AmoebaReferenceVdwForce::setSigmaCombiningRule(const std::string& sigmaCombiningRule) {
_sigmaCombiningRule = sigmaCombiningRule;
// convert to upper case and set combining function
std::transform( _sigmaCombiningRule.begin(), _sigmaCombiningRule.end(), _sigmaCombiningRule.begin(), (int(*)(int)) std::toupper);
if( _sigmaCombiningRule == "GEOMETRIC" ){
std::transform(_sigmaCombiningRule.begin(), _sigmaCombiningRule.end(), _sigmaCombiningRule.begin(), (int(*)(int)) std::toupper);
if (_sigmaCombiningRule == "GEOMETRIC") {
_combineSigmas = &AmoebaReferenceVdwForce::geometricSigmaCombiningRule;
} else if( _sigmaCombiningRule == "CUBIC-MEAN" ){
} else if (_sigmaCombiningRule == "CUBIC-MEAN") {
_combineSigmas = &AmoebaReferenceVdwForce::cubicMeanSigmaCombiningRule;
} else {
_combineSigmas = &AmoebaReferenceVdwForce::arithmeticSigmaCombiningRule;
}
}
std::string AmoebaReferenceVdwForce::getSigmaCombiningRule( void ) const {
std::string AmoebaReferenceVdwForce::getSigmaCombiningRule() const {
return _sigmaCombiningRule;
}
RealOpenMM AmoebaReferenceVdwForce::arithmeticSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const {
RealOpenMM AmoebaReferenceVdwForce::arithmeticSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const {
return (sigmaI + sigmaJ);
}
RealOpenMM AmoebaReferenceVdwForce::geometricSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const {
RealOpenMM AmoebaReferenceVdwForce::geometricSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const {
return 2.0*SQRT(sigmaI*sigmaJ);
}
RealOpenMM AmoebaReferenceVdwForce::cubicMeanSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const {
RealOpenMM AmoebaReferenceVdwForce::cubicMeanSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const {
const RealOpenMM zero = 0.0;
......@@ -120,48 +120,48 @@ RealOpenMM AmoebaReferenceVdwForce::cubicMeanSigmaCombiningRule( RealOpenMM sigm
return sigmaI != zero && sigmaJ != 0.0 ? 2.0*(sigmaI2*sigmaI + sigmaJ2*sigmaJ)/(sigmaI2 + sigmaJ2) : zero;
}
void AmoebaReferenceVdwForce::setEpsilonCombiningRule( const std::string& epsilonCombiningRule ){
void AmoebaReferenceVdwForce::setEpsilonCombiningRule(const std::string& epsilonCombiningRule) {
_epsilonCombiningRule = epsilonCombiningRule;
std::transform( _epsilonCombiningRule.begin(), _epsilonCombiningRule.end(), _epsilonCombiningRule.begin(), (int(*)(int)) std::toupper);
std::transform(_epsilonCombiningRule.begin(), _epsilonCombiningRule.end(), _epsilonCombiningRule.begin(), (int(*)(int)) std::toupper);
// convert to upper case and set combining function
if( _epsilonCombiningRule == "ARITHMETIC" ){
if (_epsilonCombiningRule == "ARITHMETIC") {
_combineEpsilons = &AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule;
} else if( _epsilonCombiningRule == "HARMONIC" ){
} else if (_epsilonCombiningRule == "HARMONIC") {
_combineEpsilons = &AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule;
} else if( _epsilonCombiningRule == "HHG" ){
} else if (_epsilonCombiningRule == "HHG") {
_combineEpsilons = &AmoebaReferenceVdwForce::hhgEpsilonCombiningRule;
} else {
_combineEpsilons = &AmoebaReferenceVdwForce::geometricEpsilonCombiningRule;
}
}
std::string AmoebaReferenceVdwForce::getEpsilonCombiningRule( void ) const {
std::string AmoebaReferenceVdwForce::getEpsilonCombiningRule() const {
return _epsilonCombiningRule;
}
RealOpenMM AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const {
RealOpenMM AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const {
return 0.5*(epsilonI + epsilonJ);
}
RealOpenMM AmoebaReferenceVdwForce::geometricEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const {
RealOpenMM AmoebaReferenceVdwForce::geometricEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const {
return SQRT(epsilonI*epsilonJ);
}
RealOpenMM AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const {
RealOpenMM AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const {
return (epsilonI != 0.0 && epsilonJ != 0.0) ? 2.0*(epsilonI*epsilonJ)/(epsilonI + epsilonJ) : 0.0;
}
RealOpenMM AmoebaReferenceVdwForce::hhgEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const {
RealOpenMM AmoebaReferenceVdwForce::hhgEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const {
RealOpenMM denominator = SQRT(epsilonI) + SQRT(epsilonJ);
return (epsilonI != 0.0 && epsilonJ != 0.0) ? 4.0*(epsilonI*epsilonJ)/(denominator*denominator) : 0.0;
}
void AmoebaReferenceVdwForce::addReducedForce( unsigned int particleI, unsigned int particleIV,
RealOpenMM reduction, RealOpenMM sign,
Vec3& force, vector<RealVec>& forces ) const {
void AmoebaReferenceVdwForce::addReducedForce(unsigned int particleI, unsigned int particleIV,
RealOpenMM reduction, RealOpenMM sign,
Vec3& force, vector<RealVec>& forces) const {
// ---------------------------------------------------------------------------------------
......@@ -178,10 +178,10 @@ void AmoebaReferenceVdwForce::addReducedForce( unsigned int particleI, unsigned
forces[particleIV][2] += sign*force[2]*(one - reduction);
}
RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combinedSigma, RealOpenMM combinedEpsilon,
const Vec3& particleIPosition,
const Vec3& particleJPosition,
Vec3& force ) const {
RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn(RealOpenMM combinedSigma, RealOpenMM combinedEpsilon,
const Vec3& particleIPosition,
const Vec3& particleJPosition,
Vec3& force) const {
// ---------------------------------------------------------------------------------------
......@@ -221,13 +221,13 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combinedSigma,
RealOpenMM ratio = (sigma_7/rho);
RealOpenMM gtau = combinedEpsilon*tau_7*r_ij_6*(ghal+one)*ratio*ratio;
RealOpenMM energy = combinedEpsilon*tau_7*sigma_7*( (ghal+one)*sigma_7/rho - two);
RealOpenMM energy = combinedEpsilon*tau_7*sigma_7*((ghal+one)*sigma_7/rho - two);
RealOpenMM dEdR = -seven*(dtau*energy + gtau);
// tapering
if( (_nonbondedMethod == CutoffNonPeriodic || _nonbondedMethod == CutoffPeriodic) && r_ij > _taperCutoff ){
if ((_nonbondedMethod == CutoffNonPeriodic || _nonbondedMethod == CutoffPeriodic) && r_ij > _taperCutoff) {
RealOpenMM delta = r_ij - _taperCutoff;
RealOpenMM taper = 1.0 + delta*delta*delta*(_taperCoefficients[C3] + delta*(_taperCoefficients[C4] + delta*_taperCoefficients[C5]));
RealOpenMM dtaper = delta*delta*(3.0*_taperCoefficients[C3] + delta*(4.0*_taperCoefficients[C4] + delta*5.0*_taperCoefficients[C5]));
......@@ -245,35 +245,35 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combinedSigma,
}
void AmoebaReferenceVdwForce::setReducedPositions( int numParticles,
const vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& reductions,
std::vector<Vec3>& reducedPositions ) const {
void AmoebaReferenceVdwForce::setReducedPositions(int numParticles,
const vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& reductions,
std::vector<Vec3>& reducedPositions) const {
static const RealOpenMM zero = 0.0;
reducedPositions.resize(numParticles);
for( unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++ ){
if( reductions[ii] != zero ){
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++) {
if (reductions[ii] != zero) {
int reductionIndex = indexIVs[ii];
reducedPositions[ii] = Vec3( reductions[ii]*( particlePositions[ii][0] - particlePositions[reductionIndex][0] ) + particlePositions[reductionIndex][0],
reductions[ii]*( particlePositions[ii][1] - particlePositions[reductionIndex][1] ) + particlePositions[reductionIndex][1],
reductions[ii]*( particlePositions[ii][2] - particlePositions[reductionIndex][2] ) + particlePositions[reductionIndex][2] );
reducedPositions[ii] = Vec3(reductions[ii]*(particlePositions[ii][0] - particlePositions[reductionIndex][0]) + particlePositions[reductionIndex][0],
reductions[ii]*(particlePositions[ii][1] - particlePositions[reductionIndex][1]) + particlePositions[reductionIndex][1],
reductions[ii]*(particlePositions[ii][2] - particlePositions[reductionIndex][2]) + particlePositions[reductionIndex][2]);
} else {
reducedPositions[ii] = Vec3( particlePositions[ii][0], particlePositions[ii][1], particlePositions[ii][2] );
reducedPositions[ii] = Vec3(particlePositions[ii][0], particlePositions[ii][1], particlePositions[ii][2]);
}
}
}
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
const vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas,
const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const std::vector< std::set<int> >& allExclusions,
vector<RealVec>& forces ) const {
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
const vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas,
const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const std::vector< std::set<int> >& allExclusions,
vector<RealVec>& forces) const {
// ---------------------------------------------------------------------------------------
......@@ -286,7 +286,7 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
// set reduced coordinates
std::vector<Vec3> reducedPositions;
setReducedPositions( numParticles, particlePositions, indexIVs, reductions, reducedPositions );
setReducedPositions(numParticles, particlePositions, indexIVs, reductions, reducedPositions);
// loop over all particle pairs
......@@ -299,44 +299,44 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
RealOpenMM energy = zero;
std::vector<unsigned int> exclusions(numParticles, 0);
for( unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++ ){
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++) {
RealOpenMM sigmaI = sigmas[ii];
RealOpenMM epsilonI = epsilons[ii];
for( std::set<int>::const_iterator jj = allExclusions[ii].begin(); jj != allExclusions[ii].end(); jj++ ){
for (std::set<int>::const_iterator jj = allExclusions[ii].begin(); jj != allExclusions[ii].end(); jj++) {
exclusions[*jj] = 1;
}
for( unsigned int jj = ii+1; jj < static_cast<unsigned int>(numParticles); jj++ ){
if( exclusions[jj] == 0 ){
for (unsigned int jj = ii+1; jj < static_cast<unsigned int>(numParticles); jj++) {
if (exclusions[jj] == 0) {
RealOpenMM combinedSigma = (this->*_combineSigmas)(sigmaI, sigmas[jj] );
RealOpenMM combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj] );
RealOpenMM combinedSigma = (this->*_combineSigmas)(sigmaI, sigmas[jj]);
RealOpenMM combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj]);
Vec3 force;
energy += calculatePairIxn( combinedSigma, combinedEpsilon,
reducedPositions[ii], reducedPositions[jj],
force );
energy += calculatePairIxn(combinedSigma, combinedEpsilon,
reducedPositions[ii], reducedPositions[jj],
force);
if( indexIVs[ii] == ii ){
if (indexIVs[ii] == ii) {
forces[ii][0] -= force[0];
forces[ii][1] -= force[1];
forces[ii][2] -= force[2];
} else {
addReducedForce( ii, indexIVs[ii], reductions[ii], -one, force, forces );
addReducedForce(ii, indexIVs[ii], reductions[ii], -one, force, forces);
}
if( indexIVs[jj] == jj ){
if (indexIVs[jj] == jj) {
forces[jj][0] += force[0];
forces[jj][1] += force[1];
forces[jj][2] += force[2];
} else {
addReducedForce( jj, indexIVs[jj], reductions[jj], one, force, forces );
addReducedForce(jj, indexIVs[jj], reductions[jj], one, force, forces);
}
}
}
for( std::set<int>::const_iterator jj = allExclusions[ii].begin(); jj != allExclusions[ii].end(); jj++ ){
for (std::set<int>::const_iterator jj = allExclusions[ii].begin(); jj != allExclusions[ii].end(); jj++) {
exclusions[*jj] = 0;
}
}
......@@ -344,14 +344,14 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
return energy;
}
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
const vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas,
const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const NeighborList& neighborList,
vector<RealVec>& forces ) const {
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
const vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas,
const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const NeighborList& neighborList,
vector<RealVec>& forces) const {
// ---------------------------------------------------------------------------------------
......@@ -364,7 +364,7 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
// set reduced coordinates
std::vector<Vec3> reducedPositions;
setReducedPositions( numParticles, particlePositions, indexIVs, reductions, reducedPositions );
setReducedPositions(numParticles, particlePositions, indexIVs, reductions, reducedPositions);
// loop over neighbor list
// (1) calculate pair vdw ixn
......@@ -373,32 +373,32 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
// based on reduction factor
RealOpenMM energy = zero;
for( unsigned int ii = 0; ii < neighborList.size(); ii++ ){
for (unsigned int ii = 0; ii < neighborList.size(); ii++) {
OpenMM::AtomPair pair = neighborList[ii];
int siteI = pair.first;
int siteJ = pair.second;
RealOpenMM combinedSigma = (this->*_combineSigmas)( sigmas[siteI], sigmas[siteJ] );
RealOpenMM combinedEpsilon = (this->*_combineEpsilons)( epsilons[siteI], epsilons[siteJ] );
RealOpenMM combinedSigma = (this->*_combineSigmas)(sigmas[siteI], sigmas[siteJ]);
RealOpenMM combinedEpsilon = (this->*_combineEpsilons)(epsilons[siteI], epsilons[siteJ]);
Vec3 force;
energy += calculatePairIxn( combinedSigma, combinedEpsilon,
reducedPositions[siteI], reducedPositions[siteJ], force );
energy += calculatePairIxn(combinedSigma, combinedEpsilon,
reducedPositions[siteI], reducedPositions[siteJ], force);
if( indexIVs[siteI] == siteI ){
if (indexIVs[siteI] == siteI) {
forces[siteI][0] -= force[0];
forces[siteI][1] -= force[1];
forces[siteI][2] -= force[2];
} else {
addReducedForce( siteI, indexIVs[siteI], reductions[siteI], -one, force, forces );
addReducedForce(siteI, indexIVs[siteI], reductions[siteI], -one, force, forces);
}
if( indexIVs[siteJ] == siteJ ){
if (indexIVs[siteJ] == siteJ) {
forces[siteJ][0] += force[0];
forces[siteJ][1] += force[1];
forces[siteJ][2] += force[2];
} else {
addReducedForce( siteJ, indexIVs[siteJ], reductions[siteJ], one, force, forces );
addReducedForce(siteJ, indexIVs[siteJ], reductions[siteJ], one, force, forces);
}
}
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.h
View file @ 41cd79a5
......@@ -34,7 +34,7 @@
namespace OpenMM {
class AmoebaReferenceVdwForce;
typedef RealOpenMM (AmoebaReferenceVdwForce::*CombiningFunction)( RealOpenMM x, RealOpenMM y) const;
typedef RealOpenMM (AmoebaReferenceVdwForce::*CombiningFunction)(RealOpenMM x, RealOpenMM y) const;
// ---------------------------------------------------------------------------------------
......@@ -71,7 +71,7 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceVdwForce( );
AmoebaReferenceVdwForce();
/**---------------------------------------------------------------------------------------
......@@ -79,8 +79,8 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceVdwForce( const std::string& sigmaCombiningRule,
const std::string& epsilonCombiningRule );
AmoebaReferenceVdwForce(const std::string& sigmaCombiningRule,
const std::string& epsilonCombiningRule);
/**---------------------------------------------------------------------------------------
......@@ -88,7 +88,7 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferenceVdwForce( ){};
~AmoebaReferenceVdwForce() {};
/**---------------------------------------------------------------------------------------
......@@ -98,7 +98,7 @@ public:
--------------------------------------------------------------------------------------- */
NonbondedMethod getNonbondedMethod( void ) const;
NonbondedMethod getNonbondedMethod() const;
/**---------------------------------------------------------------------------------------
......@@ -108,7 +108,7 @@ public:
--------------------------------------------------------------------------------------- */
void setNonbondedMethod( NonbondedMethod nonbondedMethod );
void setNonbondedMethod(NonbondedMethod nonbondedMethod);
/**---------------------------------------------------------------------------------------
......@@ -118,7 +118,7 @@ public:
--------------------------------------------------------------------------------------- */
double getCutoff( void ) const;
double getCutoff() const;
/**---------------------------------------------------------------------------------------
......@@ -128,7 +128,7 @@ public:
--------------------------------------------------------------------------------------- */
void setCutoff( double cutoff );
void setCutoff(double cutoff);
/**---------------------------------------------------------------------------------------
......@@ -138,7 +138,7 @@ public:
--------------------------------------------------------------------------------------- */
void setSigmaCombiningRule( const std::string& sigmaCombiningRule );
void setSigmaCombiningRule(const std::string& sigmaCombiningRule);
/**---------------------------------------------------------------------------------------
......@@ -148,7 +148,7 @@ public:
--------------------------------------------------------------------------------------- */
std::string getSigmaCombiningRule( void ) const;
std::string getSigmaCombiningRule() const;
/**---------------------------------------------------------------------------------------
......@@ -158,7 +158,7 @@ public:
--------------------------------------------------------------------------------------- */
void setEpsilonCombiningRule( const std::string& epsilonCombiningRule );
void setEpsilonCombiningRule(const std::string& epsilonCombiningRule);
/**---------------------------------------------------------------------------------------
......@@ -168,7 +168,7 @@ public:
--------------------------------------------------------------------------------------- */
std::string getEpsilonCombiningRule( void ) const;
std::string getEpsilonCombiningRule() const;
/**---------------------------------------------------------------------------------------
......@@ -197,12 +197,12 @@ public:
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numParticles, const std::vector<OpenMM::RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const std::vector< std::set<int> >& vdwExclusions,
std::vector<OpenMM::RealVec>& forces ) const;
RealOpenMM calculateForceAndEnergy(int numParticles, const std::vector<OpenMM::RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const std::vector< std::set<int> >& vdwExclusions,
std::vector<OpenMM::RealVec>& forces) const;
/**---------------------------------------------------------------------------------------
......@@ -221,12 +221,12 @@ public:
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numParticles, const std::vector<OpenMM::RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const NeighborList& neighborList,
std::vector<OpenMM::RealVec>& forces ) const;
RealOpenMM calculateForceAndEnergy(int numParticles, const std::vector<OpenMM::RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const NeighborList& neighborList,
std::vector<OpenMM::RealVec>& forces) const;
private:
......@@ -245,15 +245,15 @@ private:
RealOpenMM _taperCoefficients[3];
RealVec _periodicBoxVectors[3];
CombiningFunction _combineSigmas;
RealOpenMM arithmeticSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const;
RealOpenMM geometricSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const;
RealOpenMM cubicMeanSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const;
RealOpenMM arithmeticSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const;
RealOpenMM geometricSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const;
RealOpenMM cubicMeanSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const;
CombiningFunction _combineEpsilons;
RealOpenMM arithmeticEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const;
RealOpenMM geometricEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const;
RealOpenMM harmonicEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const;
RealOpenMM hhgEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const;
RealOpenMM arithmeticEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const;
RealOpenMM geometricEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const;
RealOpenMM harmonicEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ) const;
RealOpenMM hhgEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ) const;
/**---------------------------------------------------------------------------------------
......@@ -272,9 +272,9 @@ private:
--------------------------------------------------------------------------------------- */
void setReducedPositions( int numParticles, const std::vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<RealOpenMM>& reductions,
std::vector<Vec3>& reducedPositions ) const;
void setReducedPositions(int numParticles, const std::vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<RealOpenMM>& reductions,
std::vector<Vec3>& reducedPositions) const;
/**---------------------------------------------------------------------------------------
......@@ -289,9 +289,9 @@ private:
--------------------------------------------------------------------------------------- */
void addReducedForce( unsigned int particleI, unsigned int particleIV,
RealOpenMM reduction, RealOpenMM sign,
Vec3& force, std::vector<OpenMM::RealVec>& forces ) const;
void addReducedForce(unsigned int particleI, unsigned int particleIV,
RealOpenMM reduction, RealOpenMM sign,
Vec3& force, std::vector<OpenMM::RealVec>& forces) const;
/**---------------------------------------------------------------------------------------
......@@ -301,7 +301,7 @@ private:
--------------------------------------------------------------------------------------- */
void setTaperCoefficients( double cutoff );
void setTaperCoefficients(double cutoff);
/**---------------------------------------------------------------------------------------
......@@ -317,9 +317,9 @@ private:
--------------------------------------------------------------------------------------- */
RealOpenMM calculatePairIxn( RealOpenMM combindedSigma, RealOpenMM combindedEpsilon,
const Vec3& particleIPosition, const Vec3& particleJPosition,
Vec3& force ) const;
RealOpenMM calculatePairIxn(RealOpenMM combindedSigma, RealOpenMM combindedEpsilon,
const Vec3& particleIPosition, const Vec3& particleJPosition,
Vec3& force) const;
};
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceWcaDispersionForce.cpp
View file @ 41cd79a5
......@@ -28,17 +28,17 @@
using std::vector;
using namespace OpenMM;
AmoebaReferenceWcaDispersionForce::AmoebaReferenceWcaDispersionForce( RealOpenMM epso, RealOpenMM epsh, RealOpenMM rmino, RealOpenMM rminh,
RealOpenMM awater, RealOpenMM shctd, RealOpenMM dispoff, RealOpenMM slevy ) :
AmoebaReferenceWcaDispersionForce::AmoebaReferenceWcaDispersionForce(RealOpenMM epso, RealOpenMM epsh, RealOpenMM rmino, RealOpenMM rminh,
RealOpenMM awater, RealOpenMM shctd, RealOpenMM dispoff, RealOpenMM slevy) :
_epso(epso), _epsh(epsh), _rmino(rmino), _rminh(rminh), _awater(awater), _shctd(shctd), _dispoff(dispoff), _slevy(slevy) {
}
RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiusI, RealOpenMM radiusK,
const RealVec& particleIPosition,
const RealVec& particleJPosition,
const RealOpenMM* const intermediateValues,
Vec3& force ) const {
RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn(RealOpenMM radiusI, RealOpenMM radiusK,
const RealVec& particleIPosition,
const RealVec& particleJPosition,
const RealOpenMM* const intermediateValues,
Vec3& force) const {
// ---------------------------------------------------------------------------------------
......@@ -62,7 +62,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
RealOpenMM zr = particleIPosition[2] - particleJPosition[2];
RealOpenMM r2 = xr*xr + yr*yr + zr*zr;
RealOpenMM r = SQRT( r2 );
RealOpenMM r = SQRT(r2);
RealOpenMM r3 = r2*r;
RealOpenMM sK = radiusK*_shctd;
......@@ -84,7 +84,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
RealOpenMM sum = zero;
RealOpenMM de = zero;
if( radiusI < (r + sK) ){
if (radiusI < (r + sK)) {
RealOpenMM rmax = (radiusI > (r - sK)) ? radiusI : (r - sK);
......@@ -93,7 +93,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
RealOpenMM lik3 = lik2*lik;
RealOpenMM lik4 = lik2*lik2;
if( lik < rmixo ){
if (lik < rmixo) {
RealOpenMM uik = (r + sK) < rmixo ? (r + sK) : rmixo;
RealOpenMM uik2 = uik*uik;
......@@ -103,7 +103,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
RealOpenMM term = four*PI/(fortyEight*r)* (three*(lik4-uik4) - eight*r*(lik3-uik3) + six*(r2-sK2)*(lik2-uik2));
RealOpenMM dl;
if( radiusI > (r - sK) ){
if (radiusI > (r - sK)) {
dl = -lik2 + two*(r2 + sK2);
dl *= lik2;
} else {
......@@ -112,7 +112,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
}
RealOpenMM du;
if( (r+sK) > rmixo ){
if ((r+sK) > rmixo) {
du = -uik2 + two*(r2 + sK2);
du *= -uik2;
} else {
......@@ -123,7 +123,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
sum += -emixo*term;
}
if( lik < rmixh ){
if (lik < rmixh) {
RealOpenMM uik = (r + sK) < rmixh ? (r + sK) : rmixh;
RealOpenMM uik2 = uik*uik;
......@@ -131,7 +131,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
RealOpenMM uik4 = uik2*uik2;
RealOpenMM term = four*PI / (fortyEight*r)*(three*(lik4-uik4) - eight*r*(lik3-uik3) + six*(r2-sK2)*(lik2-uik2));
RealOpenMM dl;
if( radiusI > (r-sK) ){
if (radiusI > (r-sK)) {
dl = -lik2 + two*(r2 + sK2);
dl *= lik2;
} else {
......@@ -140,7 +140,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
}
RealOpenMM du;
if (r+sK > rmixh){
if (r+sK > rmixh) {
du = -uik2 + two*(r2 + sK2);
du *= -uik2;
} else {
......@@ -162,7 +162,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
RealOpenMM uik12 = uik6 * uik6;
RealOpenMM uik13 = uik10 * uik3;
if( uik > rmixo ){
if (uik > rmixo) {
RealOpenMM lik = rmax > rmixo ? rmax : rmixo;
RealOpenMM lik2 = lik * lik;
......@@ -177,7 +177,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
RealOpenMM term = four*PI/(120.0*r*lik5*uik5)*(15.0*uik*lik*r*(uik4-lik4) - ten*uik2*lik2*(uik3-lik3) + six*(sK2-r2)*(uik5-lik5));
RealOpenMM dl;
if( radiusI > (r-sK) || rmax < rmixo ){
if (radiusI > (r-sK) || rmax < rmixo) {
dl = -five*lik2 + three*(r2 + sK2);
dl /= -lik5;
} else {
......@@ -189,7 +189,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
RealOpenMM idisp = -two*ao*term;
de -= two*ao*PI*(dl + du)/(15.0*r2);
term = four*PI/(2640.0*r*lik12*uik12) * (120.0*uik*lik*r*(uik11-lik11) - 66.0*uik2*lik2*(uik10-lik10) + 55.0*(sK2-r2)*(uik12-lik12));
if( radiusI > (r-sK) || rmax < rmixo ){
if (radiusI > (r-sK) || rmax < rmixo) {
dl = -six*lik2 + five*r2 + five*sK2;
dl /= -lik12;
} else {
......@@ -201,7 +201,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
de += ao*rmixo7*PI*(dl + du)/(60.0*r2);
sum += ao*rmixo7*term + idisp;
}
if (uik > rmixh){
if (uik > rmixh) {
lik = rmax > rmixh ? rmax : rmixh;
lik2 = lik * lik;
lik3 = lik2 * lik;
......@@ -215,7 +215,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
RealOpenMM term = four*PI / (120.0*r*lik5*uik5)*(15.0*uik*lik*r*(uik4-lik4) - ten*uik2*lik2*(uik3-lik3) + six*(sK2-r2)*(uik5-lik5));
RealOpenMM dl;
if( radiusI > (r-sK) || rmax < rmixh ){
if (radiusI > (r-sK) || rmax < rmixh) {
dl = -five*lik2 + three*(r2 + sK2);
dl /= -lik5;
} else {
......@@ -227,7 +227,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
RealOpenMM idisp = -four*ah*term;
de = de - four*ah*PI*(dl + du)/(15.0*r2);
term = four*PI / (2640.0*r*lik12*uik12)* (120.0*uik*lik*r*(uik11-lik11)- 66.0*uik2*lik2*(uik10-lik10)+ 55.0*(sK2-r2)*(uik12-lik12));
if( radiusI > (r-sK) || rmax < rmixh){
if (radiusI > (r-sK) || rmax < rmixh) {
dl = -six*lik2 + five*r2 + five*sK2;
dl = -dl / lik12;
} else {;
......@@ -254,12 +254,12 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn( RealOpenMM radiu
}
RealOpenMM AmoebaReferenceWcaDispersionForce::calculateForceAndEnergy( int numParticles,
const vector<RealVec>& particlePositions,
const std::vector<RealOpenMM>& radii,
const std::vector<RealOpenMM>& epsilons,
RealOpenMM totalMaximumDispersionEnergy,
vector<RealVec>& forces ) const {
RealOpenMM AmoebaReferenceWcaDispersionForce::calculateForceAndEnergy(int numParticles,
const vector<RealVec>& particlePositions,
const std::vector<RealOpenMM>& radii,
const std::vector<RealOpenMM>& epsilons,
RealOpenMM totalMaximumDispersionEnergy,
vector<RealVec>& forces) const {
// ---------------------------------------------------------------------------------------
......@@ -284,7 +284,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculateForceAndEnergy( int numPa
RealOpenMM intermediateValues[LastIntermediateValueIndex];
for( unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++ ){
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++) {
RealOpenMM epsi = epsilons[ii];
RealOpenMM rmini = radii[ii];
......@@ -296,7 +296,7 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculateForceAndEnergy( int numPa
RealOpenMM rminI2 = rmini*rmini;
RealOpenMM rminI3 = rminI2*rmini;
RealOpenMM rmixo = two*(rmino3 + rminI3 ) / (rmino2 + rminI2);
RealOpenMM rmixo = two*(rmino3 + rminI3) / (rmino2 + rminI2);
intermediateValues[RMIXO] = rmixo;
RealOpenMM rmixo7 = rmixo*rmixo*rmixo;
......@@ -319,14 +319,14 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculateForceAndEnergy( int numPa
intermediateValues[AH] = emixh*rmixh7;
for( unsigned int jj = 0; jj < static_cast<unsigned int>(numParticles); jj++ ){
for (unsigned int jj = 0; jj < static_cast<unsigned int>(numParticles); jj++) {
if( ii == jj )continue;
if (ii == jj)continue;
Vec3 force;
energy += calculatePairIxn( rmini, radii[jj],
particlePositions[ii], particlePositions[jj],
intermediateValues, force );
energy += calculatePairIxn(rmini, radii[jj],
particlePositions[ii], particlePositions[jj],
intermediateValues, force);
forces[ii][0] += force[0];
forces[ii][1] += force[1];
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceWcaDispersionForce.h
View file @ 41cd79a5
......@@ -51,8 +51,8 @@ public:
--------------------------------------------------------------------------------------- */
AmoebaReferenceWcaDispersionForce( RealOpenMM epso, RealOpenMM epsh, RealOpenMM rmino, RealOpenMM rminh,
RealOpenMM awater, RealOpenMM shctd, RealOpenMM dispoff, RealOpenMM slevy );
AmoebaReferenceWcaDispersionForce(RealOpenMM epso, RealOpenMM epsh, RealOpenMM rmino, RealOpenMM rminh,
RealOpenMM awater, RealOpenMM shctd, RealOpenMM dispoff, RealOpenMM slevy);
/**---------------------------------------------------------------------------------------
......@@ -60,7 +60,7 @@ public:
--------------------------------------------------------------------------------------- */
~AmoebaReferenceWcaDispersionForce( ){};
~AmoebaReferenceWcaDispersionForce() {};
/**---------------------------------------------------------------------------------------
......@@ -77,10 +77,10 @@ public:
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numParticles, const std::vector<OpenMM::RealVec>& particlePositions,
const std::vector<RealOpenMM>& radii,
const std::vector<RealOpenMM>& epsilons,
RealOpenMM totalMaximumDispersionEnergy, std::vector<OpenMM::RealVec>& forces ) const;
RealOpenMM calculateForceAndEnergy(int numParticles, const std::vector<OpenMM::RealVec>& particlePositions,
const std::vector<RealOpenMM>& radii,
const std::vector<RealOpenMM>& epsilons,
RealOpenMM totalMaximumDispersionEnergy, std::vector<OpenMM::RealVec>& forces) const;
private:
RealOpenMM _epso;
......@@ -109,10 +109,10 @@ private:
--------------------------------------------------------------------------------------- */
RealOpenMM calculatePairIxn( RealOpenMM radiusI, RealOpenMM radiusJ,
const OpenMM::RealVec& particleIPosition, const OpenMM::RealVec& particleJPosition,
const RealOpenMM* const intermediateValues,
Vec3& force ) const;
RealOpenMM calculatePairIxn(RealOpenMM radiusI, RealOpenMM radiusJ,
const OpenMM::RealVec& particleIPosition, const OpenMM::RealVec& particleJPosition,
const RealOpenMM* const intermediateValues,
Vec3& force) const;
};
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaAngleForce.cpp
View file @ 41cd79a5
......@@ -66,7 +66,7 @@ const double TOL = 1e-5;
--------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){
static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
......@@ -75,24 +75,24 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return;
}
static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, double quadraticK, double cubicK,
double quarticK, double penticK, double sexticK,
double* dEdR, double* energyTerm, FILE* log ) {
static void getPrefactorsGivenAngleCosine(double cosine, double idealAngle, double quadraticK, double cubicK,
double quarticK, double penticK, double sexticK,
double* dEdR, double* energyTerm, FILE* log) {
double angle;
if( cosine >= 1.0 ){
if (cosine >= 1.0) {
angle = 0.0f;
}
else if( cosine <= -1.0 ){
else if (cosine <= -1.0) {
angle = RADIAN*PI_M;
}
else {
angle = RADIAN*acos(cosine);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "getPrefactorsGivenAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealAngle );
(void) fflush( log );
if (log) {
(void) fprintf(log, "getPrefactorsGivenAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealAngle);
(void) fflush(log);
}
#endif
......@@ -103,11 +103,11 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou
// deltaIdeal = r - r_0
*dEdR = ( 2.0 +
3.0*cubicK* deltaIdeal +
4.0*quarticK*deltaIdeal2 +
5.0*penticK* deltaIdeal3 +
6.0*sexticK* deltaIdeal4 );
*dEdR = (2.0 +
3.0*cubicK* deltaIdeal +
4.0*quarticK*deltaIdeal2 +
5.0*penticK* deltaIdeal3 +
6.0*sexticK* deltaIdeal4 );
*dEdR *= RADIAN*quadraticK*deltaIdeal;
......@@ -122,29 +122,29 @@ static void getPrefactorsGivenAngleCosine( double cosine, double idealAngle, dou
}
static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions, AmoebaAngleForce& AmoebaAngleForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) {
std::vector<Vec3>& forces, double* energy, FILE* log) {
int particle1, particle2, particle3;
double idealAngle;
double quadraticK;
AmoebaAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, idealAngle, quadraticK );
AmoebaAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, idealAngle, quadraticK);
double cubicK = AmoebaAngleForce.getAmoebaGlobalAngleCubic();
double quarticK = AmoebaAngleForce.getAmoebaGlobalAngleQuartic();
double penticK = AmoebaAngleForce.getAmoebaGlobalAnglePentic();
double sexticK = AmoebaAngleForce.getAmoebaGlobalAngleSextic();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaAngleForce: bond %d [%d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, idealAngle, quadraticK, cubicK, quarticK, penticK, sexticK );
(void) fflush( log );
if (log) {
(void) fprintf(log, "computeAmoebaAngleForce: bond %d [%d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, idealAngle, quadraticK, cubicK, quarticK, penticK, sexticK);
(void) fflush(log);
}
#endif
double deltaR[2][3];
double r2_0 = 0.0;
double r2_1 = 0.0;
for( int ii = 0; ii < 3; ii++ ){
for (int ii = 0; ii < 3; ii++) {
deltaR[0][ii] = positions[particle1][ii] - positions[particle2][ii];
r2_0 += deltaR[0][ii]*deltaR[0][ii];
......@@ -155,33 +155,33 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions
}
double pVector[3];
crossProductVector3( deltaR[0], deltaR[1], pVector );
double rp = sqrt( pVector[0]*pVector[0] + pVector[1]*pVector[1] + pVector[2]*pVector[2] );
if( rp < 1.0e-06 ){
crossProductVector3(deltaR[0], deltaR[1], pVector);
double rp = sqrt(pVector[0]*pVector[0] + pVector[1]*pVector[1] + pVector[2]*pVector[2]);
if (rp < 1.0e-06) {
rp = 1.0e-06;
}
double dot = deltaR[0][0]*deltaR[1][0] + deltaR[0][1]*deltaR[1][1] + deltaR[0][2]*deltaR[1][2];
double cosine = dot/sqrt(r2_0*r2_1);
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "dot=%10.3e r2_0=%10.3e r2_1=%10.3e\n", dot, r2_0, r2_1 );
(void) fflush( log );
if (log) {
(void) fprintf(log, "dot=%10.3e r2_0=%10.3e r2_1=%10.3e\n", dot, r2_0, r2_1);
(void) fflush(log);
}
#endif
double dEdR;
double energyTerm;
getPrefactorsGivenAngleCosine( cosine, idealAngle, quadraticK, cubicK,
quarticK, penticK, sexticK, &dEdR, &energyTerm, log );
getPrefactorsGivenAngleCosine(cosine, idealAngle, quadraticK, cubicK,
quarticK, penticK, sexticK, &dEdR, &energyTerm, log);
double termA = -dEdR/(r2_0*rp);
double termC = dEdR/(r2_1*rp);
double deltaCrossP[3][3];
crossProductVector3( deltaR[0], pVector, deltaCrossP[0] );
crossProductVector3( deltaR[1], pVector, deltaCrossP[2] );
for( int ii = 0; ii < 3; ii++ ){
crossProductVector3(deltaR[0], pVector, deltaCrossP[0]);
crossProductVector3(deltaR[1], pVector, deltaCrossP[2]);
for (int ii = 0; ii < 3; ii++) {
deltaCrossP[0][ii] *= termA;
deltaCrossP[2][ii] *= termC;
deltaCrossP[1][ii] = -1.0*(deltaCrossP[0][ii] + deltaCrossP[2][ii]);
......@@ -202,71 +202,71 @@ static void computeAmoebaAngleForce(int bondIndex, std::vector<Vec3>& positions
*energy += energyTerm;
}
static void computeAmoebaAngleForces( Context& context, AmoebaAngleForce& AmoebaAngleForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) {
static void computeAmoebaAngleForces(Context& context, AmoebaAngleForce& AmoebaAngleForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log) {
// get positions and zero forces
State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() );
expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){
for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
}
// calculates forces/energy
*expectedEnergy = 0.0;
for( int ii = 0; ii < AmoebaAngleForce.getNumAngles(); ii++ ){
computeAmoebaAngleForce(ii, positions, AmoebaAngleForce, expectedForces, expectedEnergy, log );
for (int ii = 0; ii < AmoebaAngleForce.getNumAngles(); ii++) {
computeAmoebaAngleForce(ii, positions, AmoebaAngleForce, expectedForces, expectedEnergy, log);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaAngleForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaAngleForces: expected energy=%14.7e\n", *expectedEnergy);
for (unsigned int ii = 0; ii < positions.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
return;
}
void compareWithExpectedForceAndEnergy( Context& context, AmoebaAngleForce& AmoebaAngleForce,
double tolerance, const std::string& idString, FILE* log) {
void compareWithExpectedForceAndEnergy(Context& context, AmoebaAngleForce& AmoebaAngleForce,
double tolerance, const std::string& idString, FILE* log) {
std::vector<Vec3> expectedForces;
double expectedEnergy;
computeAmoebaAngleForces( context, AmoebaAngleForce, expectedForces, &expectedEnergy, log );
computeAmoebaAngleForces(context, AmoebaAngleForce, expectedForces, &expectedEnergy, log);
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy());
for (unsigned int ii = 0; ii < forces.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
}
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
}
void testOneAngle( FILE* log ) {
void testOneAngle(FILE* log) {
System system;
int numberOfParticles = 3;
for( int ii = 0; ii < numberOfParticles; ii++ ){
for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0);
}
......@@ -290,7 +290,7 @@ void testOneAngle( FILE* log ) {
system.addForce(amoebaAngleForce);
ASSERT(!amoebaAngleForce->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, Platform::getPlatformByName( "Reference"));
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles);
......@@ -299,7 +299,7 @@ void testOneAngle( FILE* log ) {
positions[2] = Vec3(0, 0, 1);
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log );
compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle", log);
// Try changing the angle parameters and make sure it's still correct.
......@@ -307,29 +307,29 @@ void testOneAngle( FILE* log ) {
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log );
compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle", log);
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaAngleForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaAngleForce, TOL, "testOneAngle", log );
compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle", log);
}
int main( int numberOfArguments, char* argv[] ) {
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestCudaAmoebaAngleForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories();
//FILE* log = fopen( "AmoebaAngleForce.log", "w" );;
//FILE* log = fopen("AmoebaAngleForce.log", "w");;
FILE* log = NULL;
//FILE* log = stderr;
testOneAngle( log );
testOneAngle(log);
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
if (log && log != stderr)
(void) fclose(log);
#endif
}
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaBondForce.cpp
View file @ 41cd79a5
......@@ -49,22 +49,22 @@ extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();
const double TOL = 1e-5;
static void computeAmoebaBondForce(int bondIndex, std::vector<Vec3>& positions, AmoebaBondForce& AmoebaBondForce,
std::vector<Vec3>& forces, double* energy ) {
std::vector<Vec3>& forces, double* energy) {
int particle1, particle2;
double bondLength;
double quadraticK;
double cubicK = AmoebaBondForce.getAmoebaGlobalBondCubic();
double quarticK = AmoebaBondForce.getAmoebaGlobalBondQuartic();
AmoebaBondForce.getBondParameters(bondIndex, particle1, particle2, bondLength, quadraticK );
AmoebaBondForce.getBondParameters(bondIndex, particle1, particle2, bondLength, quadraticK);
double deltaR[3];
double r2 = 0.0;
for( int ii = 0; ii < 3; ii++ ){
for (int ii = 0; ii < 3; ii++) {
deltaR[ii] = positions[particle2][ii] - positions[particle1][ii];
r2 += deltaR[ii]*deltaR[ii];
}
double r = sqrt( r2 );
double r = sqrt(r2);
double bondDelta = (r - bondLength);
double bondDelta2 = bondDelta*bondDelta;
......@@ -84,64 +84,64 @@ static void computeAmoebaBondForce(int bondIndex, std::vector<Vec3>& positions,
}
static void computeAmoebaBondForces( Context& context, AmoebaBondForce& AmoebaBondForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) {
static void computeAmoebaBondForces(Context& context, AmoebaBondForce& AmoebaBondForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log) {
// get positions and zero forces
State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() );
expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){
for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
}
// calculates forces/energy
*expectedEnergy = 0.0;
for( int ii = 0; ii < AmoebaBondForce.getNumBonds(); ii++ ){
computeAmoebaBondForce(ii, positions, AmoebaBondForce, expectedForces, expectedEnergy );
for (int ii = 0; ii < AmoebaBondForce.getNumBonds(); ii++) {
computeAmoebaBondForce(ii, positions, AmoebaBondForce, expectedForces, expectedEnergy);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaBondForces: expected energy=%15.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%15.7e %15.7e %15.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaBondForces: expected energy=%15.7e\n", *expectedEnergy);
for (unsigned int ii = 0; ii < positions.size(); ii++) {
(void) fprintf(log, "%6u [%15.7e %15.7e %15.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
return;
}
void compareWithExpectedForceAndEnergy( Context& context, AmoebaBondForce& AmoebaBondForce, double tolerance, const std::string& idString, FILE* log) {
void compareWithExpectedForceAndEnergy(Context& context, AmoebaBondForce& AmoebaBondForce, double tolerance, const std::string& idString, FILE* log) {
std::vector<Vec3> expectedForces;
double expectedEnergy;
computeAmoebaBondForces( context, AmoebaBondForce, expectedForces, &expectedEnergy, NULL );
computeAmoebaBondForces(context, AmoebaBondForce, expectedForces, &expectedEnergy, NULL);
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaBondForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaBondForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy());
for (unsigned int ii = 0; ii < forces.size(); ii++) {
(void) fprintf(log, "%6u [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
}
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
}
void testOneBond( FILE* log ) {
void testOneBond(FILE* log) {
System system;
......@@ -156,22 +156,22 @@ void testOneBond( FILE* log ) {
double quadraticK = 1.0;
double cubicK = 2.0;
double quarticicK = 3.0;
amoebaBondForce->setAmoebaGlobalBondCubic( cubicK );
amoebaBondForce->setAmoebaGlobalBondQuartic( quarticicK );
amoebaBondForce->setAmoebaGlobalBondCubic(cubicK);
amoebaBondForce->setAmoebaGlobalBondQuartic(quarticicK);
amoebaBondForce->addBond(0, 1, bondLength, quadraticK);
system.addForce(amoebaBondForce);
Context context(system, integrator, Platform::getPlatformByName( "Reference"));
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(2);
positions[0] = Vec3(0, 1, 0);
positions[1] = Vec3(0, 0, 0);
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testOneBond", log );
compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testOneBond", log);
}
void testTwoBond( FILE* log ) {
void testTwoBond(FILE* log) {
System system;
......@@ -187,15 +187,15 @@ void testTwoBond( FILE* log ) {
double quadraticK = 1.0;
double cubicK = 2.0;
double quarticicK = 3.0;
amoebaBondForce->setAmoebaGlobalBondCubic( cubicK );
amoebaBondForce->setAmoebaGlobalBondQuartic( quarticicK );
amoebaBondForce->setAmoebaGlobalBondCubic(cubicK);
amoebaBondForce->setAmoebaGlobalBondQuartic(quarticicK);
amoebaBondForce->addBond(0, 1, bondLength, quadraticK);
amoebaBondForce->addBond(1, 2, bondLength, quadraticK);
system.addForce(amoebaBondForce);
ASSERT(!amoebaBondForce->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, Platform::getPlatformByName( "Reference"));
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(3);
positions[0] = Vec3(0, 1, 0);
......@@ -203,7 +203,7 @@ void testTwoBond( FILE* log ) {
positions[2] = Vec3(1, 0, 1);
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log );
compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond", log);
// Try changing the bond parameters and make sure it's still correct.
......@@ -212,17 +212,17 @@ void testTwoBond( FILE* log ) {
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log );
compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond", log);
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaBondForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaBondForce, TOL, "testTwoBond", log );
compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond", log);
}
int main( int numberOfArguments, char* argv[] ) {
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestReferenceAmoebaBondForce running test..." << std::endl;
......@@ -230,11 +230,11 @@ int main( int numberOfArguments, char* argv[] ) {
FILE* log = NULL;
//FILE* log = stderr;
//testOneBond( log );
testTwoBond( log );
//testOneBond(log);
testTwoBond(log);
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
if (log && log != stderr)
(void) fclose(log);
#endif
}
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaGeneralizedKirkwoodForce.cpp
View file @ 41cd79a5
......@@ -65,10 +65,10 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
int numberOfParticles = 8;
 
amoebaMultipoleForce->setNonbondedMethod( AmoebaMultipoleForce::NoCutoff );
amoebaMultipoleForce->setPolarizationType( polarizationType );
amoebaMultipoleForce->setMutualInducedTargetEpsilon( 1.0e-06 );
amoebaMultipoleForce->setMutualInducedMaxIterations( 500 );
amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::NoCutoff);
amoebaMultipoleForce->setPolarizationType(polarizationType);
amoebaMultipoleForce->setMutualInducedTargetEpsilon(1.0e-06);
amoebaMultipoleForce->setMutualInducedMaxIterations(500);
 
std::vector<double> nitrogenMolecularDipole(3);
std::vector<double> nitrogenMolecularQuadrupole(9);
......@@ -89,8 +89,8 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
 
// first N
 
system.addParticle( 1.4007000e+01 );
amoebaMultipoleForce->addMultipole( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 1, 2, 3, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03 );
system.addParticle(1.4007000e+01);
amoebaMultipoleForce->addMultipole( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 1, 2, 3, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03);
 
// 3 H attached to first N
 
......@@ -110,167 +110,167 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
hydrogenMolecularQuadrupole[7] = 0.0000000e+00;
hydrogenMolecularQuadrupole[8] = 2.4549167e-06;
 
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 2, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 2, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
system.addParticle(1.0080000e+00);
system.addParticle(1.0080000e+00);
system.addParticle(1.0080000e+00);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 2, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 2, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
 
// second N
 
system.addParticle( 1.4007000e+01 );
amoebaMultipoleForce->addMultipole( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 5, 6, 7, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03 );
system.addParticle( 1.4007000e+01);
amoebaMultipoleForce->addMultipole( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 5, 6, 7, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03);
 
// 3 H attached to second N
 
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 6, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 6, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
system.addParticle( 1.0080000e+00);
system.addParticle( 1.0080000e+00);
system.addParticle( 1.0080000e+00);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 6, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
amoebaMultipoleForce->addMultipole( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 6, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04);
 
// covalent maps
 
std::vector< int > covalentMap;
covalentMap.resize(0);
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.push_back(1);
covalentMap.push_back(2);
covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap(0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 0 );
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.push_back(0);
covalentMap.push_back(1);
covalentMap.push_back(2);
covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap(0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 0 );
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.push_back(0);
amoebaMultipoleForce->setCovalentMap(1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.push_back(2);
covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap(1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 0 );
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.push_back(0);
covalentMap.push_back(1);
covalentMap.push_back(2);
covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap(1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 0 );
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.push_back(0);
amoebaMultipoleForce->setCovalentMap(2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 1 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.push_back(1);
covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap(2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 0 );
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.push_back(0);
covalentMap.push_back(1);
covalentMap.push_back(2);
covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap(2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 0 );
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.push_back(0);
amoebaMultipoleForce->setCovalentMap(3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.push_back(1);
covalentMap.push_back(2);
amoebaMultipoleForce->setCovalentMap(3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 0 );
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.push_back(0);
covalentMap.push_back(1);
covalentMap.push_back(2);
covalentMap.push_back(3);
amoebaMultipoleForce->setCovalentMap(3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.push_back(5);
covalentMap.push_back(6);
covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap(4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 4 );
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.push_back(4);
covalentMap.push_back(5);
covalentMap.push_back(6);
covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap(4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 4 );
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.push_back(4);
amoebaMultipoleForce->setCovalentMap(5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.push_back(6);
covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap(5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 4 );
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.push_back(4);
covalentMap.push_back(5);
covalentMap.push_back(6);
covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap(5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 4 );
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.push_back(4);
amoebaMultipoleForce->setCovalentMap(6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 5 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.push_back(5);
covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap(6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 4 );
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.push_back(4);
covalentMap.push_back(5);
covalentMap.push_back(6);
covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap(6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 4 );
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.push_back(4);
amoebaMultipoleForce->setCovalentMap(7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.push_back(5);
covalentMap.push_back(6);
amoebaMultipoleForce->setCovalentMap(7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap);
 
covalentMap.resize(0);
covalentMap.push_back( 4 );
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.push_back(4);
covalentMap.push_back(5);
covalentMap.push_back(6);
covalentMap.push_back(7);
amoebaMultipoleForce->setCovalentMap(7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap);
 
system.addForce(amoebaMultipoleForce);
 
// GK force
 
amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01 );
amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00 );
amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm( includeCavityTerm );
amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01);
amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00);
amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm(includeCavityTerm);
 
// addParticle: charge, radius, scalingFactor
 
for( unsigned int ii = 0; ii < 2; ii++ ){
amoebaGeneralizedKirkwoodForce->addParticle( -5.7960000e-01, 1.5965000e-01, 6.9000000e-01 );
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 );
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 );
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 );
for (unsigned int ii = 0; ii < 2; ii++) {
amoebaGeneralizedKirkwoodForce->addParticle( -5.7960000e-01, 1.5965000e-01, 6.9000000e-01);
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01);
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01);
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01);
}
system.addForce(amoebaGeneralizedKirkwoodForce);
}
......@@ -278,14 +278,14 @@ static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce
static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>& forces, double& energy, FILE* log) {
std::vector<Vec3> positions(context.getSystem().getNumParticles());
 
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03 );
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02 );
positions[2] = Vec3( 2.0843610e-01, 8.0953200e-02, 3.7462200e-02 );
positions[3] = Vec3( 1.7280780e-01, 2.0730000e-04, -9.8741700e-02 );
positions[4] = Vec3( -1.6743680e-01, 1.5900000e-05, -6.6149000e-03 );
positions[5] = Vec3( -2.0428260e-01, 8.1071500e-02, 4.1343900e-02 );
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02 );
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02 );
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03);
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02);
positions[2] = Vec3( 2.0843610e-01, 8.0953200e-02, 3.7462200e-02);
positions[3] = Vec3( 1.7280780e-01, 2.0730000e-04, -9.8741700e-02);
positions[4] = Vec3( -1.6743680e-01, 1.5900000e-05, -6.6149000e-03);
positions[5] = Vec3( -2.0428260e-01, 8.1071500e-02, 4.1343900e-02);
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02);
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02);
 
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
......@@ -295,8 +295,8 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>&
 
// setup for villin
 
static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::PolarizationType polarizationType,
int includeCavityTerm, std::vector<Vec3>& forces, double& energy, FILE* log ){
static void setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::PolarizationType polarizationType,
int includeCavityTerm, std::vector<Vec3>& forces, double& energy, FILE* log) {
 
// beginning of Multipole setup
 
......@@ -305,13 +305,13 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
int numberOfParticles = 596;
 
amoebaMultipoleForce->setNonbondedMethod( AmoebaMultipoleForce::NoCutoff );
amoebaMultipoleForce->setPolarizationType( polarizationType );
amoebaMultipoleForce->setMutualInducedTargetEpsilon( 1.0e-06 );
amoebaMultipoleForce->setMutualInducedMaxIterations( 500 );
amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::NoCutoff);
amoebaMultipoleForce->setPolarizationType(polarizationType);
amoebaMultipoleForce->setMutualInducedTargetEpsilon(1.0e-06);
amoebaMultipoleForce->setMutualInducedMaxIterations(500);
 
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){
system.addParticle( 1.0 );
for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0);
}
 
static const double multipoleData[] = {
......@@ -929,7 +929,7 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
std::vector<double> quadrupole(9);
unsigned int entriesPerParticle = 21;
const double* data = multipoleData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){
for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
 
dipole[0] = data[dipoleIndex + 0];
dipole[1] = data[dipoleIndex + 1];
......@@ -945,9 +945,9 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
quadrupole[7] = data[quadrupoleIndex + 7];
quadrupole[8] = data[quadrupoleIndex + 8];
 
amoebaMultipoleForce->addMultipole( data[chargeIndex], dipole, quadrupole, static_cast<int>(data[axisTypeIndex]),
amoebaMultipoleForce->addMultipole(data[chargeIndex], dipole, quadrupole, static_cast<int>(data[axisTypeIndex]),
static_cast<int>(data[multipoleAtomZIndex]), static_cast<int>(data[multipoleAtomXIndex]), static_cast<int>(data[multipoleAtomYIndex]),
data[tholeIndex], data[dampingFactorIndex], data[polarityIndex] );
data[tholeIndex], data[dampingFactorIndex], data[polarityIndex]);
data += entriesPerParticle;
}
 
......@@ -5727,15 +5727,15 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
unsigned int covalentMapDataSize = sizeof(covalentMapData)/sizeof(int);
 
unsigned int covalentIndex = 0;
while( covalentIndex < covalentMapDataSize ){
while (covalentIndex < covalentMapDataSize) {
int particleIndex = covalentMapData[covalentIndex++];
int typeIndex = covalentMapData[covalentIndex++];
int entries = covalentMapData[covalentIndex++];
std::vector< int > covalentMap(entries);
for( unsigned int ii = 0; ii < entries; ii++ ){
for (unsigned int ii = 0; ii < entries; ii++) {
covalentMap[ii] = covalentMapData[covalentIndex++];
}
amoebaMultipoleForce->setCovalentMap( particleIndex, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(typeIndex), covalentMap );
amoebaMultipoleForce->setCovalentMap(particleIndex, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(typeIndex), covalentMap);
}
system.addForce(amoebaMultipoleForce);
 
......@@ -5744,9 +5744,9 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
// GK force
 
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01 );
amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00 );
amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm( includeCavityTerm );
amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01);
amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00);
amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm(includeCavityTerm);
 
// addParticle: charge, radius, scalingFactor
 
......@@ -6350,8 +6350,8 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
};
 
const double* gkData = generalizedKirkwoodData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){
amoebaGeneralizedKirkwoodForce->addParticle( gkData[0], gkData[1], gkData[2] );
for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
amoebaGeneralizedKirkwoodForce->addParticle(gkData[0], gkData[1], gkData[2]);
gkData += 3;
}
system.addForce(amoebaGeneralizedKirkwoodForce);
......@@ -6960,15 +6960,15 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
};
 
const double* positionsDataPtr = positionsData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){
positions[ii] = Vec3( positionsDataPtr[0], positionsDataPtr[1], positionsDataPtr[2] );
for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
positions[ii] = Vec3(positionsDataPtr[0], positionsDataPtr[1], positionsDataPtr[2]);
positionsDataPtr += 3;
}
 
std::string platformName;
platformName = "Reference";
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName( platformName ) );
Context context(system, integrator, Platform::getPlatformByName(platformName));
 
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
......@@ -6978,117 +6978,117 @@ static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Polari
 
// compare forces and energies
 
static void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy,
const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) {
static void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy,
const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log) {
 
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
if( log ){
if (log) {
double conversion = 1.0/4.184;
double energyAbsDiff = fabs( expectedEnergy - energy );
double energyRelDiff = 2.0*energyAbsDiff/( fabs( expectedEnergy ) + fabs( energy ) + 1.0e-08 );
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e absDiff=%15.7e relDiff=%15.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy,
conversion*energyAbsDiff, conversion*energyRelDiff );
if( conversion != 1.0 )conversion *= -0.1;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2] );
double norm = sqrt( forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2] );
double absDiff = fabs( norm - expectedNorm );
double relDiff = 2.0*absDiff/(fabs( norm ) + fabs( expectedNorm ) + 1.0e-08);
(void) fprintf( log, "%6u %15.7e %15.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*absDiff, conversion*relDiff,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2], conversion*expectedNorm, conversion*norm );
double energyAbsDiff = fabs(expectedEnergy - energy);
double energyRelDiff = 2.0*energyAbsDiff/(fabs(expectedEnergy) + fabs(energy) + 1.0e-08);
(void) fprintf(log, "%s: expected energy=%14.7e %14.7e absDiff=%15.7e relDiff=%15.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy,
conversion*energyAbsDiff, conversion*energyRelDiff);
if (conversion != 1.0)conversion *= -0.1;
for (unsigned int ii = 0; ii < forces.size(); ii++) {
double expectedNorm = sqrt(expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2]);
double norm = sqrt(forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2]);
double absDiff = fabs(norm - expectedNorm);
double relDiff = 2.0*absDiff/(fabs(norm) + fabs(expectedNorm) + 1.0e-08);
(void) fprintf(log, "%6u %15.7e %15.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*absDiff, conversion*relDiff,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2], conversion*expectedNorm, conversion*norm);
}
(void) fflush( log );
(void) fflush(log);
conversion = 1.0;
(void) fprintf( log, "\n%s: expected energy=%14.7e %14.7e no conversion\n", testName.c_str(), conversion*expectedEnergy, conversion*energy );
if( conversion != 1.0 )conversion = -1.0;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2] );
(void) fprintf(log, "\n%s: expected energy=%14.7e %14.7e no conversion\n", testName.c_str(), conversion*expectedEnergy, conversion*energy);
if (conversion != 1.0)conversion = -1.0;
for (unsigned int ii = 0; ii < forces.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
 
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC_MOD( expectedForces[ii], forces[ii], tolerance, testName );
for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC_MOD(expectedForces[ii], forces[ii], tolerance, testName);
}
ASSERT_EQUAL_TOL_MOD( expectedEnergy, energy, tolerance, testName );
ASSERT_EQUAL_TOL_MOD(expectedEnergy, energy, tolerance, testName);
}
 
// compare relative differences in force norms and energies
 
static void compareForceNormsEnergy( std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) {
static void compareForceNormsEnergy(std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log) {
 
 
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
if( log ){
if (log) {
double conversion = 1.0/4.184;
double energyAbsDiff = fabs( expectedEnergy - energy );
double energyRelDiff = 2.0*energyAbsDiff/( fabs( expectedEnergy ) + fabs( energy ) + 1.0e-08 );
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e absDiff=%15.7e relDiff=%15.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy,
conversion*energyAbsDiff, conversion*energyRelDiff );
if( conversion != 1.0 )conversion *= -0.1;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2] );
double norm = sqrt( forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2] );
double absDiff = fabs( (norm - expectedNorm) );
double relDiff = 2.0*absDiff/(fabs( norm ) + fabs( expectedNorm ) + 1.0e-08);
(void) fprintf( log, "%6u %15.7e %15.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e] %15.7e %15.7e\n", ii,
double energyAbsDiff = fabs(expectedEnergy - energy);
double energyRelDiff = 2.0*energyAbsDiff/(fabs(expectedEnergy) + fabs(energy) + 1.0e-08);
(void) fprintf(log, "%s: expected energy=%14.7e %14.7e absDiff=%15.7e relDiff=%15.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy,
conversion*energyAbsDiff, conversion*energyRelDiff);
if (conversion != 1.0)conversion *= -0.1;
for (unsigned int ii = 0; ii < forces.size(); ii++) {
double expectedNorm = sqrt(expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2]);
double norm = sqrt(forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2]);
double absDiff = fabs((norm - expectedNorm));
double relDiff = 2.0*absDiff/(fabs(norm) + fabs(expectedNorm) + 1.0e-08);
(void) fprintf(log, "%6u %15.7e %15.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e] %15.7e %15.7e\n", ii,
fabs(conversion)*absDiff, relDiff,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2],
fabs(conversion)*expectedNorm, fabs(conversion)*norm );
fabs(conversion)*expectedNorm, fabs(conversion)*norm);
}
(void) fflush( log );
(void) fflush(log);
conversion = 1.0;
(void) fprintf( log, "\n%s: expected energy=%14.7e %14.7e no conversion\n", testName.c_str(), conversion*expectedEnergy, conversion*energy );
if( conversion != 1.0 )conversion = -1.0;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
(void) fprintf(log, "\n%s: expected energy=%14.7e %14.7e no conversion\n", testName.c_str(), conversion*expectedEnergy, conversion*energy);
if (conversion != 1.0)conversion = -1.0;
for (unsigned int ii = 0; ii < forces.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2] );
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
 
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2] );
for (unsigned int ii = 0; ii < forces.size(); ii++) {
double expectedNorm = sqrt(expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2]);
 
double norm = sqrt( forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2] );
double absDiff = fabs( norm - expectedNorm );
double relDiff = 2.0*absDiff/(fabs( norm ) + fabs( expectedNorm ) + 1.0e-08);
double norm = sqrt(forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2]);
double absDiff = fabs(norm - expectedNorm);
double relDiff = 2.0*absDiff/(fabs(norm) + fabs(expectedNorm) + 1.0e-08);
 
if( relDiff > tolerance && absDiff > 0.001 ){
if (relDiff > tolerance && absDiff > 0.001) {
std::stringstream details;
details << testName << "Relative difference in norms " << relDiff << " larger than allowed tolerance at particle=" << ii;
details << ": norms=" << norm << " expected norm=" << expectedNorm;
throwException(__FILE__, __LINE__, details.str());
}
}
double energyAbsDiff = fabs( expectedEnergy - energy );
double energyRelDiff = 2.0*energyAbsDiff/( fabs( expectedEnergy ) + fabs( energy ) + 1.0e-08 );
if( energyRelDiff > tolerance ){
double energyAbsDiff = fabs(expectedEnergy - energy);
double energyRelDiff = 2.0*energyAbsDiff/(fabs(expectedEnergy) + fabs(energy) + 1.0e-08);
if (energyRelDiff > tolerance) {
std::stringstream details;
details << testName << "Relative difference in energies " << energyRelDiff << " larger than allowed tolerance.";
details << "Energies=" << energy << " expected energy=" << expectedEnergy;
......@@ -7098,7 +7098,7 @@ static void compareForceNormsEnergy( std::string& testName, double expectedEnerg
 
// test GK direct polarization for system comprised of two ammonia molecules
 
static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) {
static void testGeneralizedKirkwoodAmmoniaDirectPolarization(FILE* log) {
 
std::string testName = "testGeneralizedKirkwoodAmmoniaDirectPolarization";
 
......@@ -7111,27 +7111,27 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) {
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log );
getForcesEnergyMultipoleAmmonia(context, forces, energy, log);
std::vector<Vec3> expectedForces(numberOfParticles);
 
double expectedEnergy = -7.6636680e+01;
 
expectedForces[0] = Vec3( -6.9252994e+02, -8.9085133e+00, 9.6489739e+01 );
expectedForces[1] = Vec3( 1.5593797e+02, -6.0331931e+01, 1.5104507e+01 );
expectedForces[2] = Vec3( 1.5870088e+02, 6.1702809e+01, 6.7708985e+00 );
expectedForces[3] = Vec3( 1.4089885e+02, 7.5870617e+00, -1.1362294e+02 );
expectedForces[4] = Vec3( -1.8916205e+02, 2.1465549e-01, -4.3433152e+02 );
expectedForces[5] = Vec3( 1.0208290e+01, 6.2676753e+01, 1.4987953e+02 );
expectedForces[6] = Vec3( 4.0621859e+02, 1.8962203e-01, 1.3021956e+02 );
expectedForces[7] = Vec3( 9.7274235e+00, -6.3130458e+01, 1.4949024e+02 );
expectedForces[0] = Vec3( -6.9252994e+02, -8.9085133e+00, 9.6489739e+01);
expectedForces[1] = Vec3( 1.5593797e+02, -6.0331931e+01, 1.5104507e+01);
expectedForces[2] = Vec3( 1.5870088e+02, 6.1702809e+01, 6.7708985e+00);
expectedForces[3] = Vec3( 1.4089885e+02, 7.5870617e+00, -1.1362294e+02);
expectedForces[4] = Vec3( -1.8916205e+02, 2.1465549e-01, -4.3433152e+02);
expectedForces[5] = Vec3( 1.0208290e+01, 6.2676753e+01, 1.4987953e+02);
expectedForces[6] = Vec3( 4.0621859e+02, 1.8962203e-01, 1.3021956e+02);
expectedForces[7] = Vec3( 9.7274235e+00, -6.3130458e+01, 1.4949024e+02);
 
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log);
}
 
// test GK mutual polarization for system comprised of two ammonia molecules
 
static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) {
static void testGeneralizedKirkwoodAmmoniaMutualPolarization(FILE* log) {
 
std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarization";
 
......@@ -7144,28 +7144,28 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) {
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Mutual, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log );
getForcesEnergyMultipoleAmmonia(context, forces, energy, log);
std::vector<Vec3> expectedForces(numberOfParticles);
 
double expectedEnergy = -7.8018875e+01;
 
expectedForces[0] = Vec3( -7.6820301e+02, -1.0102760e+01, 1.0094389e+02 );
expectedForces[1] = Vec3( 1.7037307e+02, -7.5621857e+01, 2.3320365e+01 );
expectedForces[2] = Vec3( 1.7353828e+02, 7.7199741e+01, 1.3965379e+01 );
expectedForces[3] = Vec3( 1.5045244e+02, 8.5784569e+00, -1.3377619e+02 );
expectedForces[4] = Vec3( -2.1811615e+02, -1.6818022e-01, -4.6103163e+02 );
expectedForces[5] = Vec3( 6.2091942e+00, 7.6748687e+01, 1.5883463e+02 );
expectedForces[6] = Vec3( 4.8035662e+02, 4.9704902e-01, 1.3948083e+02 );
expectedForces[7] = Vec3( 5.3895456e+00, -7.7131137e+01, 1.5826273e+02 );
expectedForces[0] = Vec3( -7.6820301e+02, -1.0102760e+01, 1.0094389e+02);
expectedForces[1] = Vec3( 1.7037307e+02, -7.5621857e+01, 2.3320365e+01);
expectedForces[2] = Vec3( 1.7353828e+02, 7.7199741e+01, 1.3965379e+01);
expectedForces[3] = Vec3( 1.5045244e+02, 8.5784569e+00, -1.3377619e+02);
expectedForces[4] = Vec3( -2.1811615e+02, -1.6818022e-01, -4.6103163e+02);
expectedForces[5] = Vec3( 6.2091942e+00, 7.6748687e+01, 1.5883463e+02);
expectedForces[6] = Vec3( 4.8035662e+02, 4.9704902e-01, 1.3948083e+02);
expectedForces[7] = Vec3( 5.3895456e+00, -7.7131137e+01, 1.5826273e+02);
 
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log);
}
 
// test GK mutual polarization for system comprised of two ammonia molecules
// including cavity term
 
static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE* log ) {
static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(FILE* log) {
 
std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm";
 
......@@ -7180,22 +7180,22 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE
ASSERT(!system.usesPeriodicBoundaryConditions());
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
getForcesEnergyMultipoleAmmonia(context, forces, energy, log );
getForcesEnergyMultipoleAmmonia(context, forces, energy, log);
std::vector<Vec3> expectedForces(numberOfParticles);
 
double expectedEnergy = -6.0434582e+01;
 
expectedForces[0] = Vec3( -7.8323218e+02, -1.0097644e+01, 1.0256890e+02 );
expectedForces[1] = Vec3( 1.7078480e+02, -7.1896701e+01, 2.0840172e+01 );
expectedForces[2] = Vec3( 1.7394089e+02, 7.3488594e+01, 1.1484648e+01 );
expectedForces[3] = Vec3( 1.5169364e+02, 8.5611299e+00, -1.2968050e+02 );
expectedForces[4] = Vec3( -2.1669693e+02, -1.5926823e-01, -4.6636274e+02 );
expectedForces[5] = Vec3( 8.7397444e+00, 7.3330990e+01, 1.6016898e+02 );
expectedForces[6] = Vec3( 4.8684950e+02, 4.8937161e-01, 1.4137061e+02 );
expectedForces[7] = Vec3( 7.9205382e+00, -7.3716473e+01, 1.5960993e+02 );
expectedForces[0] = Vec3( -7.8323218e+02, -1.0097644e+01, 1.0256890e+02);
expectedForces[1] = Vec3( 1.7078480e+02, -7.1896701e+01, 2.0840172e+01);
expectedForces[2] = Vec3( 1.7394089e+02, 7.3488594e+01, 1.1484648e+01);
expectedForces[3] = Vec3( 1.5169364e+02, 8.5611299e+00, -1.2968050e+02);
expectedForces[4] = Vec3( -2.1669693e+02, -1.5926823e-01, -4.6636274e+02);
expectedForces[5] = Vec3( 8.7397444e+00, 7.3330990e+01, 1.6016898e+02);
expectedForces[6] = Vec3( 4.8684950e+02, 4.8937161e-01, 1.4137061e+02);
expectedForces[7] = Vec3( 7.9205382e+00, -7.3716473e+01, 1.5960993e+02);
 
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log);
// Try changing the particle parameters and make sure it's still correct.
......@@ -7225,7 +7225,7 @@ static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE
 
// test GK direct polarization for villin system
 
static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) {
static void testGeneralizedKirkwoodVillinDirectPolarization(FILE* log) {
 
std::string testName = "testGeneralizedKirkwoodVillinDirectPolarization";
 
......@@ -7233,7 +7233,7 @@ static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) {
std::vector<Vec3> forces;
double energy;
 
setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Direct, 0, forces, energy, log );
setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Direct, 0, forces, energy, log);
std::vector<Vec3> expectedForces(numberOfParticles);
 
double expectedEnergy = -8.4281157e+03;
......@@ -7838,18 +7838,18 @@ static void testGeneralizedKirkwoodVillinDirectPolarization( FILE* log ) {
};
 
const double* forceDataPtr = forceData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){
expectedForces[ii] = Vec3( forceDataPtr[0], forceDataPtr[1], forceDataPtr[2] );
for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
expectedForces[ii] = Vec3(forceDataPtr[0], forceDataPtr[1], forceDataPtr[2]);
forceDataPtr += 3;
}
 
double tolerance = 1.0e-05;
compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
compareForceNormsEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log);
}
 
// test GK mutual polarization for villin system
 
static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) {
static void testGeneralizedKirkwoodVillinMutualPolarization(FILE* log) {
 
std::string testName = "testGeneralizedKirkwoodVillinMutualPolarization";
 
......@@ -7857,7 +7857,7 @@ static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) {
std::vector<Vec3> forces;
double energy;
 
setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::Mutual, 0, forces, energy, log );
setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Mutual, 0, forces, energy, log);
std::vector<Vec3> expectedForces(numberOfParticles);
 
double expectedEnergy = -8.6477811e+03;
......@@ -8462,16 +8462,16 @@ static void testGeneralizedKirkwoodVillinMutualPolarization( FILE* log ) {
};
 
const double* forceDataPtr = forceData;
for( unsigned int ii = 0; ii < numberOfParticles; ii++ ){
expectedForces[ii] = Vec3( forceDataPtr[0], forceDataPtr[1], forceDataPtr[2] );
for (unsigned int ii = 0; ii < numberOfParticles; ii++) {
expectedForces[ii] = Vec3(forceDataPtr[0], forceDataPtr[1], forceDataPtr[2]);
forceDataPtr += 3;
}
 
double tolerance = 1.0e-05;
compareForceNormsEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
compareForceNormsEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log);
}
 
int main( int numberOfArguments, char* argv[] ) {
int main(int numberOfArguments, char* argv[]) {
 
try {
std::cout << "TestReferenceAmoebaGeneralizedKirkwoodForce running test..." << std::endl;
......@@ -8482,11 +8482,11 @@ int main( int numberOfArguments, char* argv[] ) {
// test direct and mutual polarization cases and
// mutual polarization w/ the cavity term
 
testGeneralizedKirkwoodAmmoniaMutualPolarization( log );
testGeneralizedKirkwoodAmmoniaDirectPolarization( log );
testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( log );
testGeneralizedKirkwoodVillinDirectPolarization( log );
testGeneralizedKirkwoodVillinMutualPolarization( log );
testGeneralizedKirkwoodAmmoniaMutualPolarization(log);
testGeneralizedKirkwoodAmmoniaDirectPolarization(log);
testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm(log);
testGeneralizedKirkwoodVillinDirectPolarization(log);
testGeneralizedKirkwoodVillinMutualPolarization(log);
 
}
catch(const std::exception& e) {
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaInPlaneAngleForce.cpp
View file @ 41cd79a5
......@@ -64,7 +64,7 @@ const double TOL = 1e-5;
--------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){
static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
......@@ -73,28 +73,28 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return;
}
static double dotVector3( double* vectorX, double* vectorY ){
static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
}
static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInPlaneAngle, double quadraticK, double cubicK,
double quarticK, double penticK, double sexticK,
double* dEdR, double* energyTerm, FILE* log ) {
static void getPrefactorsGivenInPlaneAngleCosine(double cosine, double idealInPlaneAngle, double quadraticK, double cubicK,
double quarticK, double penticK, double sexticK,
double* dEdR, double* energyTerm, FILE* log) {
double angle;
if( cosine >= 1.0 ){
if (cosine >= 1.0) {
angle = 0.0f;
}
else if( cosine <= -1.0 ){
else if (cosine <= -1.0) {
angle = RADIAN*PI_M;
}
else {
angle = RADIAN*acos(cosine);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "getPrefactorsGivenInPlaneAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealInPlaneAngle );
(void) fflush( log );
if (log) {
(void) fprintf(log, "getPrefactorsGivenInPlaneAngleCosine: cosine=%10.3e angle=%10.3e ideal=%10.3e\n", cosine, angle, idealInPlaneAngle);
(void) fflush(log);
}
#endif
......@@ -105,11 +105,11 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP
// deltaIdeal = r - r_0
*dEdR = ( 2.0 +
*dEdR = (2.0 +
3.0*cubicK* deltaIdeal +
4.0*quarticK*deltaIdeal2 +
5.0*penticK* deltaIdeal3 +
6.0*sexticK* deltaIdeal4 );
6.0*sexticK* deltaIdeal4 );
*dEdR *= RADIAN*quadraticK*deltaIdeal;
......@@ -124,22 +124,22 @@ static void getPrefactorsGivenInPlaneAngleCosine( double cosine, double idealInP
}
static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& positions, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) {
std::vector<Vec3>& forces, double* energy, FILE* log) {
int particle1, particle2, particle3, particle4;
double idealInPlaneAngle;
double quadraticK;
AmoebaInPlaneAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK );
AmoebaInPlaneAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK);
double cubicK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleCubic();
double quarticK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleQuartic();
double penticK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAnglePentic();
double sexticK = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleSextic();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaInPlaneAngleForce: bond %d [%d %d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK, cubicK, quarticK, penticK, sexticK );
(void) fflush( log );
if (log) {
(void) fprintf(log, "computeAmoebaInPlaneAngleForce: bond %d [%d %d %d %d] ang=%10.3f k=%10.3f [%10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK, cubicK, quarticK, penticK, sexticK);
(void) fflush(log);
}
#endif
......@@ -162,82 +162,82 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
// APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex
double deltaR[LastDeltaAtomIndex][3];
for( int ii = 0; ii < 3; ii++ ){
for (int ii = 0; ii < 3; ii++) {
deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii];
deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii];
deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii];
}
crossProductVector3( deltaR[AD], deltaR[CD], deltaR[T] );
crossProductVector3(deltaR[AD], deltaR[CD], deltaR[T]);
double rT2 = dotVector3( deltaR[T], deltaR[T] );
double delta = dotVector3( deltaR[T], deltaR[BD] );
double rT2 = dotVector3(deltaR[T], deltaR[T]);
double delta = dotVector3(deltaR[T], deltaR[BD]);
delta *= -1.0/rT2;
for( int ii = 0; ii < 3; ii++ ){
for (int ii = 0; ii < 3; ii++) {
deltaR[P][ii] = positions[particle2][ii] + deltaR[T][ii]*delta;
deltaR[AP][ii] = positions[particle1][ii] - deltaR[P][ii];
deltaR[CP][ii] = positions[particle3][ii] - deltaR[P][ii];
}
double rAp2 = dotVector3( deltaR[AP], deltaR[AP] );
double rCp2 = dotVector3( deltaR[CP], deltaR[CP] );
if( rAp2 <= 0.0 && rCp2 <= 0.0 ){
double rAp2 = dotVector3(deltaR[AP], deltaR[AP]);
double rCp2 = dotVector3(deltaR[CP], deltaR[CP]);
if (rAp2 <= 0.0 && rCp2 <= 0.0) {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaInPlaneAngleForce: rAp2 or rCp2 <= 0.0\n" );
(void) fflush( log );
if (log) {
(void) fprintf(log, "computeAmoebaInPlaneAngleForce: rAp2 or rCp2 <= 0.0\n");
(void) fflush(log);
}
#endif
return;
}
crossProductVector3( deltaR[CP], deltaR[AP], deltaR[M] );
crossProductVector3(deltaR[CP], deltaR[AP], deltaR[M]);
double rm = dotVector3( deltaR[M], deltaR[M] );
rm = sqrt( rm );
if( rm < 0.000001 ){
double rm = dotVector3(deltaR[M], deltaR[M]);
rm = sqrt(rm);
if (rm < 0.000001) {
rm = 0.000001;
}
double dot = dotVector3( deltaR[AP], deltaR[CP] );
double cosine = dot/sqrt( rAp2*rCp2 );
double dot = dotVector3(deltaR[AP], deltaR[CP]);
double cosine = dot/sqrt(rAp2*rCp2);
double dEdR;
double energyTerm;
getPrefactorsGivenInPlaneAngleCosine( cosine, idealInPlaneAngle, quadraticK, cubicK,
quarticK, penticK, sexticK, &dEdR, &energyTerm, log );
getPrefactorsGivenInPlaneAngleCosine(cosine, idealInPlaneAngle, quadraticK, cubicK,
quarticK, penticK, sexticK, &dEdR, &energyTerm, log);
double termA = -dEdR/(rAp2*rm);
double termC = dEdR/(rCp2*rm);
crossProductVector3( deltaR[AP], deltaR[M], deltaR[APxM] );
crossProductVector3( deltaR[CP], deltaR[M], deltaR[CPxM] );
crossProductVector3(deltaR[AP], deltaR[M], deltaR[APxM]);
crossProductVector3(deltaR[CP], deltaR[M], deltaR[CPxM]);
// forces will be gathered here
enum { dA, dB, dC, dD, LastDIndex };
double forceTerm[LastDIndex][3];
for( int ii = 0; ii < 3; ii++ ){
for (int ii = 0; ii < 3; ii++) {
forceTerm[dA][ii] = deltaR[APxM][ii]*termA;
forceTerm[dC][ii] = deltaR[CPxM][ii]*termC;
forceTerm[dB][ii] = -1.0*( forceTerm[dA][ii] + forceTerm[dC][ii] );
forceTerm[dB][ii] = -1.0*(forceTerm[dA][ii] + forceTerm[dC][ii]);
}
double pTrT2 = dotVector3( forceTerm[dB], deltaR[T] );
double pTrT2 = dotVector3(forceTerm[dB], deltaR[T]);
pTrT2 /= rT2;
crossProductVector3( deltaR[CD], forceTerm[dB], deltaR[CDxdB] );
crossProductVector3( forceTerm[dB], deltaR[AD], deltaR[dBxAD] );
crossProductVector3(deltaR[CD], forceTerm[dB], deltaR[CDxdB]);
crossProductVector3(forceTerm[dB], deltaR[AD], deltaR[dBxAD]);
if( fabs( pTrT2 ) > 1.0e-08 ){
if (fabs(pTrT2) > 1.0e-08) {
double delta2 = delta*2.0;
crossProductVector3( deltaR[BD], deltaR[CD], deltaR[BDxCD] );
crossProductVector3( deltaR[T], deltaR[CD], deltaR[TxCD] );
crossProductVector3( deltaR[AD], deltaR[BD], deltaR[ADxBD] );
crossProductVector3( deltaR[AD], deltaR[T], deltaR[ADxT] );
for( int ii = 0; ii < 3; ii++ ){
crossProductVector3(deltaR[BD], deltaR[CD], deltaR[BDxCD]);
crossProductVector3(deltaR[T], deltaR[CD], deltaR[TxCD] );
crossProductVector3(deltaR[AD], deltaR[BD], deltaR[ADxBD]);
crossProductVector3(deltaR[AD], deltaR[T], deltaR[ADxT] );
for (int ii = 0; ii < 3; ii++) {
double term = deltaR[BDxCD][ii] + delta2*deltaR[TxCD][ii];
forceTerm[dA][ii] += delta*deltaR[CDxdB][ii] + term*pTrT2;
......@@ -245,16 +245,16 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
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] );
forceTerm[dD][ii] = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
}
}
else {
for( int ii = 0; ii < 3; ii++ ){
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[dD][ii] = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
}
}
......@@ -280,69 +280,69 @@ static void computeAmoebaInPlaneAngleForce(int bondIndex, std::vector<Vec3>& po
}
static void computeAmoebaInPlaneAngleForces( Context& context, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) {
static void computeAmoebaInPlaneAngleForces(Context& context, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log) {
// get positions and zero forces
State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() );
expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){
for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
}
// calculates forces/energy
*expectedEnergy = 0.0;
for( int ii = 0; ii < AmoebaInPlaneAngleForce.getNumAngles(); ii++ ){
computeAmoebaInPlaneAngleForce(ii, positions, AmoebaInPlaneAngleForce, expectedForces, expectedEnergy, log );
for (int ii = 0; ii < AmoebaInPlaneAngleForce.getNumAngles(); ii++) {
computeAmoebaInPlaneAngleForce(ii, positions, AmoebaInPlaneAngleForce, expectedForces, expectedEnergy, log);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaInPlaneAngleForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaInPlaneAngleForces: expected energy=%14.7e\n", *expectedEnergy);
for (unsigned int ii = 0; ii < positions.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
return;
}
void compareWithExpectedForceAndEnergy( Context& context, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce,
double tolerance, const std::string& idString, FILE* log) {
void compareWithExpectedForceAndEnergy(Context& context, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce,
double tolerance, const std::string& idString, FILE* log) {
std::vector<Vec3> expectedForces;
double expectedEnergy;
computeAmoebaInPlaneAngleForces( context, AmoebaInPlaneAngleForce, expectedForces, &expectedEnergy, log );
computeAmoebaInPlaneAngleForces(context, AmoebaInPlaneAngleForce, expectedForces, &expectedEnergy, log);
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaInPlaneAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaInPlaneAngleForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy());
for (unsigned int ii = 0; ii < forces.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
}
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
}
void testOneAngle( FILE* log ) {
void testOneAngle(FILE* log) {
System system;
int numberOfParticles = 4;
for( int ii = 0; ii < numberOfParticles; ii++ ){
for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0);
}
......@@ -366,7 +366,7 @@ void testOneAngle( FILE* log ) {
system.addForce(amoebaInPlaneAngleForce);
ASSERT(!amoebaInPlaneAngleForce->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, Platform::getPlatformByName( "Reference"));
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles);
......@@ -376,7 +376,7 @@ void testOneAngle( FILE* log ) {
positions[3] = Vec3(1, 1, 1);
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log );
compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log);
// Try changing the angle parameters and make sure it's still correct.
......@@ -384,29 +384,29 @@ void testOneAngle( FILE* log ) {
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log );
compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log);
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaInPlaneAngleForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log );
compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle", log);
}
int main( int numberOfArguments, char* argv[] ) {
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestReferenceAmoebaInPlaneAngleForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories();
FILE* log = NULL;
//FILE* log = stderr;
//FILE* log = fopen( "AmoebaInPlaneAngleForce.log", "w" );;
//FILE* log = fopen("AmoebaInPlaneAngleForce.log", "w");;
testOneAngle( NULL );
testOneAngle(NULL);
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
if (log && log != stderr)
(void) fclose(log);
#endif
}
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaMultipoleForce.cpp
View file @ 41cd79a5
This source diff could not be displayed because it is too large. You can view the blob instead.
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaOutOfPlaneBendForce.cpp
View file @ 41cd79a5
......@@ -64,7 +64,7 @@ const double TOL = 1e-3;
--------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){
static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
......@@ -73,13 +73,13 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return;
}
static double dotVector3( double* vectorX, double* vectorY ){
static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
}
static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) {
std::vector<Vec3>& forces, double* energy, FILE* log) {
double kAngleCubic = amoebaOutOfPlaneBendForce.getAmoebaGlobalOutOfPlaneBendCubic();
......@@ -89,12 +89,12 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
int particle1, particle2, particle3, particle4;
double kAngleQuadratic;
amoebaOutOfPlaneBendForce.getOutOfPlaneBendParameters(bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic );
amoebaOutOfPlaneBendForce.getOutOfPlaneBendParameters(bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic);
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForce: bond %d [%d %d %d %d] k=[%10.3e %10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic, kAngleCubic, kAngleQuartic, kAnglePentic, kAngleSextic );
(void) fflush( log );
if (log) {
(void) fprintf(log, "computeAmoebaOutOfPlaneBendForce: bond %d [%d %d %d %d] k=[%10.3e %10.3e %10.3e %10.3e %10.3e]\n",
bondIndex, particle1, particle2, particle3, particle4, kAngleQuadratic, kAngleCubic, kAngleQuartic, kAnglePentic, kAngleSextic);
(void) fflush(log);
}
#endif
......@@ -108,7 +108,7 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
// and various intermediate terms
double deltaR[LastDeltaIndex][6];
for( int ii = 0; ii < 3; ii++ ){
for (int ii = 0; ii < 3; ii++) {
deltaR[AB][ii] = positions[particle1][ii] - positions[particle2][ii];
deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii];
deltaR[DB][ii] = positions[particle4][ii] - positions[particle2][ii];
......@@ -116,36 +116,36 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii];
}
double rDB2 = dotVector3( deltaR[DB], deltaR[DB] );
double rAD2 = dotVector3( deltaR[AD], deltaR[AD] );
double rCD2 = dotVector3( deltaR[CD], deltaR[CD] );
double rDB2 = dotVector3(deltaR[DB], deltaR[DB]);
double rAD2 = dotVector3(deltaR[AD], deltaR[AD]);
double rCD2 = dotVector3(deltaR[CD], deltaR[CD]);
double tempVector[3];
crossProductVector3( deltaR[CB], deltaR[DB], tempVector );
double eE = dotVector3( deltaR[AB], tempVector );
double dot = dotVector3( deltaR[AD], deltaR[CD] );
crossProductVector3(deltaR[CB], deltaR[DB], tempVector);
double eE = dotVector3(deltaR[AB], tempVector );
double dot = dotVector3(deltaR[AD], deltaR[CD]);
double cc = rAD2*rCD2 - dot*dot;
if( rDB2 <= 0.0 || cc == 0.0 ){
if (rDB2 <= 0.0 || cc == 0.0) {
return;
}
double bkk2 = rDB2 - eE*eE/cc;
double cosine = sqrt(bkk2/rDB2);
double angle;
if( cosine >= 1.0 ){
if (cosine >= 1.0) {
angle = 0.0;
}
else if( cosine <= -1.0 ){
else if (cosine <= -1.0) {
angle = PI_M;
}
else {
angle = RADIAN*acos(cosine);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForce: bkk2=%14.7e rDB2=%14.7e cos=%14.7e dt=%14.7e]\n",
bkk2, rDB2, cosine, angle );
(void) fflush( log );
if (log) {
(void) fprintf(log, "computeAmoebaOutOfPlaneBendForce: bkk2=%14.7e rDB2=%14.7e cos=%14.7e dt=%14.7e]\n",
bkk2, rDB2, cosine, angle);
(void) fflush(log);
}
#endif
......@@ -161,23 +161,23 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
double dEdCos;
dEdCos = dEdDt/sqrt(cc*bkk2);
if( eE > 0.0 ){
if (eE > 0.0) {
dEdCos *= -1.0;
}
double term = eE/cc;
double dccd[LastAtomIndex][3];
for( int ii = 0; ii < 3; ii++ ){
for (int ii = 0; ii < 3; ii++) {
dccd[A][ii] = (deltaR[AD][ii]*rCD2 - deltaR[CD][ii]*dot)*term;
dccd[C][ii] = (deltaR[CD][ii]*rAD2 - deltaR[AD][ii]*dot)*term;
dccd[D][ii] = -1.0*(dccd[A][ii] + dccd[C][ii]);
}
double deed[LastAtomIndex][3];
crossProductVector3( deltaR[DB], deltaR[CB], deed[A] );
crossProductVector3( deltaR[AB], deltaR[DB], deed[C] );
crossProductVector3( deltaR[CB], deltaR[AB], deed[D] );
crossProductVector3(deltaR[DB], deltaR[CB], deed[A]);
crossProductVector3(deltaR[AB], deltaR[DB], deed[C]);
crossProductVector3(deltaR[CB], deltaR[AB], deed[D]);
term = eE/rDB2;
deed[D][0] += deltaR[DB][0]*term;
......@@ -189,24 +189,24 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
// forces
// calculate forces for atoms a, c, d
// the force for b is then -( a+ c + d)
// the force for b is then -(a+ c + d)
double subForce[LastAtomIndex][3];
for( int jj = 0; jj < LastAtomIndex; jj++ ){
for (int jj = 0; jj < LastAtomIndex; jj++) {
// A, C, D
for( int ii = 0; ii < 3; ii++ ){
subForce[jj][ii] = dEdCos*( dccd[jj][ii] + deed[jj][ii] );
for (int ii = 0; ii < 3; ii++) {
subForce[jj][ii] = dEdCos*(dccd[jj][ii] + deed[jj][ii]);
}
if( jj == 0 )jj++; // skip B
if (jj == 0)jj++; // skip B
// now compute B
if( jj == 3 ){
for( int ii = 0; ii < 3; ii++ ){
if (jj == 3) {
for (int ii = 0; ii < 3; ii++) {
subForce[1][ii] = -1.0*(subForce[0][ii] + subForce[2][ii] + subForce[3][ii]);
}
}
......@@ -243,70 +243,70 @@ static void computeAmoebaOutOfPlaneBendForce(int bondIndex, std::vector<Vec3>&
return;
}
static void computeAmoebaOutOfPlaneBendForces( Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) {
static void computeAmoebaOutOfPlaneBendForces(Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log) {
// get positions and zero forces
State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() );
expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){
for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
}
// calculates forces/energy
*expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaOutOfPlaneBendForce.getNumOutOfPlaneBends(); ii++ ){
computeAmoebaOutOfPlaneBendForce(ii, positions, amoebaOutOfPlaneBendForce, expectedForces, expectedEnergy, log );
for (int ii = 0; ii < amoebaOutOfPlaneBendForce.getNumOutOfPlaneBends(); ii++) {
computeAmoebaOutOfPlaneBendForce(ii, positions, amoebaOutOfPlaneBendForce, expectedForces, expectedEnergy, log);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e\n", *expectedEnergy);
for (unsigned int ii = 0; ii < positions.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
return;
}
void compareWithExpectedForceAndEnergy( Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
double tolerance, const std::string& idString, FILE* log) {
void compareWithExpectedForceAndEnergy(Context& context, AmoebaOutOfPlaneBendForce& amoebaOutOfPlaneBendForce,
double tolerance, const std::string& idString, FILE* log) {
std::vector<Vec3> expectedForces;
double expectedEnergy;
computeAmoebaOutOfPlaneBendForces( context, amoebaOutOfPlaneBendForce, expectedForces, &expectedEnergy, log );
computeAmoebaOutOfPlaneBendForces(context, amoebaOutOfPlaneBendForce, expectedForces, &expectedEnergy, log);
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaOutOfPlaneBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy());
for (unsigned int ii = 0; ii < forces.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
}
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
}
void testOneOutOfPlaneBend( FILE* log ) {
void testOneOutOfPlaneBend(FILE* log) {
System system;
int numberOfParticles = 4;
for( int ii = 0; ii < numberOfParticles; ii++ ){
for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0);
}
......@@ -314,29 +314,29 @@ void testOneOutOfPlaneBend( FILE* log ) {
AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce();
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01 );
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic( 0.5600000E-04 );
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic( -0.7000000E-06 );
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07 );
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->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend);
system.addForce(amoebaOutOfPlaneBendForce);
ASSERT(!amoebaOutOfPlaneBendForce->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, Platform::getPlatformByName( "Reference"));
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 0.262660000E+02, 0.254130000E+02, 0.284200000E+01 );
positions[1] = Vec3( 0.269130000E+02, 0.266390000E+02, 0.353100000E+01 );
positions[0] = Vec3(0.262660000E+02, 0.254130000E+02, 0.284200000E+01);
positions[1] = Vec3(0.269130000E+02, 0.266390000E+02, 0.353100000E+01);
positions[2] = Vec3( 0.278860000E+02, 0.264630000E+02, 0.426300000E+01 );
positions[3] = Vec3( 0.245568230E+02, 0.250215290E+02, 0.796852800E+01 );
positions[2] = Vec3(0.278860000E+02, 0.264630000E+02, 0.426300000E+01);
positions[3] = Vec3(0.245568230E+02, 0.250215290E+02, 0.796852800E+01);
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log );
compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log);
// Try changing the bend parameters and make sure it's still correct.
......@@ -344,21 +344,21 @@ void testOneOutOfPlaneBend( FILE* log ) {
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log );
compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log);
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaOutOfPlaneBendForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log );
compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log);
}
void testOneOutOfPlaneBend2( FILE* log, int setId ) {
void testOneOutOfPlaneBend2(FILE* log, int setId) {
System system;
int numberOfParticles = 4;
for( int ii = 0; ii < numberOfParticles; ii++ ){
for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0);
}
......@@ -366,10 +366,10 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce();
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01 );
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic( 0.5600000E-04 );
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic( -0.7000000E-06 );
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07 );
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic( 0.5600000E-04);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic( -0.7000000E-06);
amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07);
/*
285 441 442 443 444 0.328682196E-01
286 441 442 444 443 0.164493407E-01
......@@ -389,147 +389,147 @@ void testOneOutOfPlaneBend2( FILE* log, int setId ) {
*/
std::map<int,Vec3> coordinates;
coordinates[440] = Vec3( 0.893800000E+01, 0.439800000E+01, 0.343100000E+01 );
coordinates[441] = Vec3( 0.779100000E+01, 0.614600000E+01, 0.390100000E+01 );
coordinates[442] = Vec3( 0.915400000E+01, 0.683900000E+01, 0.389400000E+01 );
coordinates[443] = Vec3( 0.101770000E+02, 0.619000000E+01, 0.379900000E+01 );
coordinates[444] = Vec3( 0.921000000E+01, 0.813800000E+01, 0.398600000E+01 );
coordinates[445] = Vec3( 0.708500000E+01, 0.672900000E+01, 0.332700000E+01 );
coordinates[446] = Vec3( 0.744300000E+01, 0.605200000E+01, 0.491900000E+01 );
coordinates[447] = Vec3( 0.100820000E+02, 0.859300000E+01, 0.398200000E+01 );
coordinates[448] = Vec3( 0.838000000E+01, 0.866100000E+01, 0.406000000E+01 );
coordinates[440] = Vec3( 0.893800000E+01, 0.439800000E+01, 0.343100000E+01);
coordinates[441] = Vec3( 0.779100000E+01, 0.614600000E+01, 0.390100000E+01);
coordinates[442] = Vec3( 0.915400000E+01, 0.683900000E+01, 0.389400000E+01);
coordinates[443] = Vec3( 0.101770000E+02, 0.619000000E+01, 0.379900000E+01);
coordinates[444] = Vec3( 0.921000000E+01, 0.813800000E+01, 0.398600000E+01);
coordinates[445] = Vec3( 0.708500000E+01, 0.672900000E+01, 0.332700000E+01);
coordinates[446] = Vec3( 0.744300000E+01, 0.605200000E+01, 0.491900000E+01);
coordinates[447] = Vec3( 0.100820000E+02, 0.859300000E+01, 0.398200000E+01);
coordinates[448] = Vec3( 0.838000000E+01, 0.866100000E+01, 0.406000000E+01);
double kOutOfPlaneBend = 0.328682196E-01;
std::vector<int> particleIndices;
if( setId == 1 ){
particleIndices.push_back( 441 );
particleIndices.push_back( 442 );
particleIndices.push_back( 443 );
particleIndices.push_back( 444 );
if (setId == 1) {
particleIndices.push_back(441);
particleIndices.push_back(442);
particleIndices.push_back(443);
particleIndices.push_back(444);
kOutOfPlaneBend = 0.328682196E-01;
}
else if( setId == 2 ){
particleIndices.push_back( 441 );
particleIndices.push_back( 442 );
particleIndices.push_back( 444 );
particleIndices.push_back( 443 );
else if (setId == 2) {
particleIndices.push_back(441);
particleIndices.push_back(442);
particleIndices.push_back(444);
particleIndices.push_back(443);
kOutOfPlaneBend = 0.164493407E-01;
}
else if( setId == 3 ){
particleIndices.push_back( 443 );
particleIndices.push_back( 442 );
particleIndices.push_back( 444 );
particleIndices.push_back( 441 );
else if (setId == 3) {
particleIndices.push_back(443);
particleIndices.push_back(442);
particleIndices.push_back(444);
particleIndices.push_back(441);
kOutOfPlaneBend = 0.636650407E-02;
}
else if( setId == 4 ){
particleIndices.push_back( 442 );
particleIndices.push_back( 444 );
particleIndices.push_back( 447 );
particleIndices.push_back( 448 );
else if (setId == 4) {
particleIndices.push_back(442);
particleIndices.push_back(444);
particleIndices.push_back(447);
particleIndices.push_back(448);
kOutOfPlaneBend = 0.392956472E-02;
}
else if( setId == 5 ){
particleIndices.push_back( 442 );
particleIndices.push_back( 444 );
particleIndices.push_back( 448 );
particleIndices.push_back( 447 );
else if (setId == 5) {
particleIndices.push_back(442);
particleIndices.push_back(444);
particleIndices.push_back(448);
particleIndices.push_back(447);
kOutOfPlaneBend = 0.392956472E-02;
}
else if( setId == 6 ){
particleIndices.push_back( 447 );
particleIndices.push_back( 444 );
particleIndices.push_back( 448 );
particleIndices.push_back( 442 );
else if (setId == 6) {
particleIndices.push_back(447);
particleIndices.push_back(444);
particleIndices.push_back(448);
particleIndices.push_back(442);
kOutOfPlaneBend = 0.214755281E-01;
}
else {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "Set id %d not recognized.\n", setId );
if (log) {
(void) fprintf(log, "Set id %d not recognized.\n", setId);
}
#endif
std::stringstream buffer;
buffer << "Set id " << setId << " not recognized.";
throw OpenMMException( buffer.str() );
throw OpenMMException(buffer.str());
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "Set id %d.\n", setId );
if (log) {
(void) fprintf(log, "Set id %d.\n", setId);
}
#endif
amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend );
amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend);
system.addForce(amoebaOutOfPlaneBendForce);
Context context(system, integrator, Platform::getPlatformByName( "Reference"));
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles);
for( unsigned int ii = 0; ii < static_cast<unsigned int>(numberOfParticles); ii++ ){
if( coordinates.find( particleIndices[ii] ) == coordinates.end() ){
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numberOfParticles); ii++) {
if (coordinates.find(particleIndices[ii]) == coordinates.end()) {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "Coordinates %d not loaded.", particleIndices[ii] );
if (log) {
(void) fprintf(log, "Coordinates %d not loaded.", particleIndices[ii]);
}
#endif
std::stringstream buffer;
buffer << "Coordinates " << particleIndices[ii] << " not loaded.";
throw OpenMMException( buffer.str() );
throw OpenMMException(buffer.str());
}
positions[ii] = coordinates[particleIndices[ii]];
}
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log );
compareWithExpectedForceAndEnergy(context, *amoebaOutOfPlaneBendForce, TOL, "testOneOutOfPlaneBend", log);
static int iter = 0;
static std::map<int,Vec3> totalForces;
static double totalEnergy;
if( iter == 0 ){
totalForces[441] = Vec3( 0.0, 0.0, 0.0 );
totalForces[442] = Vec3( 0.0, 0.0, 0.0 );
totalForces[443] = Vec3( 0.0, 0.0, 0.0 );
totalForces[444] = Vec3( 0.0, 0.0, 0.0 );
totalForces[445] = Vec3( 0.0, 0.0, 0.0 );
totalForces[446] = Vec3( 0.0, 0.0, 0.0 );
totalForces[447] = Vec3( 0.0, 0.0, 0.0 );
totalForces[448] = Vec3( 0.0, 0.0, 0.0 );
totalForces[449] = Vec3( 0.0, 0.0, 0.0 );
if (iter == 0) {
totalForces[441] = Vec3( 0.0, 0.0, 0.0);
totalForces[442] = Vec3( 0.0, 0.0, 0.0);
totalForces[443] = Vec3( 0.0, 0.0, 0.0);
totalForces[444] = Vec3( 0.0, 0.0, 0.0);
totalForces[445] = Vec3( 0.0, 0.0, 0.0);
totalForces[446] = Vec3( 0.0, 0.0, 0.0);
totalForces[447] = Vec3( 0.0, 0.0, 0.0);
totalForces[448] = Vec3( 0.0, 0.0, 0.0);
totalForces[449] = Vec3( 0.0, 0.0, 0.0);
totalEnergy = 0.0;
}
iter++;
std::vector<Vec3> forces;
forces.resize( numberOfParticles );
forces.resize(numberOfParticles);
double energy;
computeAmoebaOutOfPlaneBendForce( 0, positions, *amoebaOutOfPlaneBendForce, forces, &energy, log );
computeAmoebaOutOfPlaneBendForce(0, positions, *amoebaOutOfPlaneBendForce, forces, &energy, log);
totalEnergy += energy;
for( unsigned int ii = 0; ii < static_cast<unsigned int>(numberOfParticles); ii++ ){
for( unsigned int jj = 0; jj < 3; jj++ ){
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numberOfParticles); ii++) {
for (unsigned int jj = 0; jj < 3; jj++) {
totalForces[particleIndices[ii]][jj] += forces[ii][jj];
}
}
if( iter == 6 ){
if (iter == 6) {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaOutOfPlaneBendForces: energy=%14.7e\n", totalEnergy);
for( std::map<int,Vec3>::iterator ii = totalForces.begin(); ii != totalForces.end(); ii++ ){
if (log) {
(void) fprintf(log, "computeAmoebaOutOfPlaneBendForces: energy=%14.7e\n", totalEnergy);
for (std::map<int,Vec3>::iterator ii = totalForces.begin(); ii != totalForces.end(); ii++) {
int particleIndex = ii->first;
Vec3 forces = ii->second;
(void) fprintf( log, "%6d [%14.7e %14.7e %14.7e] \n", particleIndex,
forces[0], forces[1], forces[2] );
(void) fprintf(log, "%6d [%14.7e %14.7e %14.7e] \n", particleIndex,
forces[0], forces[1], forces[2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
}
}
int main( int numberOfArguments, char* argv[] ) {
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestReferenceAmoebaOutOfPlaneBendForce running test..." << std::endl;
......@@ -537,16 +537,16 @@ int main( int numberOfArguments, char* argv[] ) {
//FILE* log = stderr;
FILE* log = NULL;
//FILE* log = fopen( "AmoebaOutOfPlaneBendForce.log", "w" );;
//FILE* log = fopen("AmoebaOutOfPlaneBendForce.log", "w");;
testOneOutOfPlaneBend( log );
//testOneOutOfPlaneBend2( log, atoi( argv[1] ) );
//for( int ii = 1; ii <= 6; ii++ ){
// testOneOutOfPlaneBend2( log, ii );
testOneOutOfPlaneBend(log);
//testOneOutOfPlaneBend2(log, atoi(argv[1]));
//for (int ii = 1; ii <= 6; ii++) {
// testOneOutOfPlaneBend2(log, ii);
//}
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
if (log && log != stderr)
(void) fclose(log);
#endif
}
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaPiTorsionForce.cpp
View file @ 41cd79a5
......@@ -63,7 +63,7 @@ const double TOL = 1e-5;
--------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){
static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
......@@ -72,23 +72,23 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return;
}
static double dotVector3( double* vectorX, double* vectorY ){
static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
}
static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& positions, AmoebaPiTorsionForce& amoebaPiTorsionForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) {
std::vector<Vec3>& forces, double* energy, FILE* log) {
int particle1, particle2, particle3, particle4, particle5, particle6;
double kTorsion;
amoebaPiTorsionForce.getPiTorsionParameters(bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion);
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForce: bond %d [%d %d %d %d %d %d] k=%10.3e\n",
bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion );
(void) fflush( log );
if (log) {
(void) fprintf(log, "computeAmoebaPiTorsionForce: bond %d [%d %d %d %d %d %d] k=%10.3e\n",
bondIndex, particle1, particle2, particle3, particle4, particle5, particle6, kTorsion);
(void) fflush(log);
}
#endif
......@@ -98,16 +98,16 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
enum { A, B, C, D, E, F, LastAtomIndex };
double d[LastAtomIndex][3];
for( int ii = 0; ii < 3; ii++ ){
for (int ii = 0; ii < 3; ii++) {
deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii];
deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii];
deltaR[EC][ii] = positions[particle5][ii] - positions[particle3][ii];
deltaR[FC][ii] = positions[particle6][ii] - positions[particle3][ii];
}
crossProductVector3( deltaR[AD], deltaR[BD], deltaR[P] );
crossProductVector3( deltaR[EC], deltaR[FC], deltaR[Q] );
for( int ii = 0; ii < 3; ii++ ){
crossProductVector3(deltaR[AD], deltaR[BD], deltaR[P]);
crossProductVector3(deltaR[EC], deltaR[FC], deltaR[Q]);
for (int ii = 0; ii < 3; ii++) {
deltaR[CP][ii] = -deltaR[P][ii];
deltaR[DC][ii] = positions[particle4][ii] - positions[particle3][ii];
deltaR[QD][ii] = deltaR[Q][ii];
......@@ -115,25 +115,25 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
deltaR[P][ii] += positions[particle3][ii];
deltaR[Q][ii] += positions[particle4][ii];
}
crossProductVector3( deltaR[CP], deltaR[DC], deltaR[T] );
crossProductVector3( deltaR[DC], deltaR[QD], deltaR[U] );
crossProductVector3( deltaR[T], deltaR[U], deltaR[TU] );
crossProductVector3(deltaR[CP], deltaR[DC], deltaR[T] );
crossProductVector3(deltaR[DC], deltaR[QD], deltaR[U] );
crossProductVector3(deltaR[T], deltaR[U], deltaR[TU]);
double rT2 = dotVector3( deltaR[T], deltaR[T] );
double rU2 = dotVector3( deltaR[U], deltaR[U] );
double rTrU = sqrt( rT2*rU2 );
if( rTrU <= 0.0 ){
double rT2 = dotVector3(deltaR[T], deltaR[T]);
double rU2 = dotVector3(deltaR[U], deltaR[U]);
double rTrU = sqrt(rT2*rU2);
if (rTrU <= 0.0) {
return;
}
double rDC = dotVector3( deltaR[DC], deltaR[DC] );
rDC = sqrt( rDC );
double rDC = dotVector3(deltaR[DC], deltaR[DC]);
rDC = sqrt(rDC);
double cosine = dotVector3( deltaR[T], deltaR[U] );
double cosine = dotVector3(deltaR[T], deltaR[U]);
cosine /= rTrU;
double sine = dotVector3( deltaR[DC], deltaR[TU] );
sine /= ( rDC*rTrU );
double sine = dotVector3(deltaR[DC], deltaR[TU]);
sine /= (rDC*rTrU);
double cosine2 = cosine*cosine - sine*sine;
double sine2 = 2.0*cosine*sine;
......@@ -143,37 +143,37 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
double dedphi = kTorsion*dphi2;
for( int ii = 0; ii < 3; ii++ ){
for (int ii = 0; ii < 3; ii++) {
deltaR[DP][ii] = positions[particle4][ii] - deltaR[P][ii];
deltaR[QC][ii] = deltaR[Q][ii] - positions[particle3][ii];
}
double factorT = dedphi/( rDC*rT2 );
double factorU = -dedphi/( rDC*rU2 );
double factorT = dedphi/(rDC*rT2);
double factorU = -dedphi/(rDC*rU2);
crossProductVector3( deltaR[T], deltaR[DC], deltaR[dT] );
crossProductVector3( deltaR[U], deltaR[DC], deltaR[dU] );
for( int ii = 0; ii < 3; ii++ ){
crossProductVector3(deltaR[T], deltaR[DC], deltaR[dT] );
crossProductVector3(deltaR[U], deltaR[DC], deltaR[dU] );
for (int ii = 0; ii < 3; ii++) {
deltaR[dT][ii] *= factorT;
deltaR[dU][ii] *= factorU;
}
crossProductVector3( deltaR[dT], deltaR[DC], deltaR[dP] );
crossProductVector3( deltaR[dU], deltaR[DC], deltaR[dQ] );
crossProductVector3(deltaR[dT], deltaR[DC], deltaR[dP] );
crossProductVector3(deltaR[dU], deltaR[DC], deltaR[dQ] );
crossProductVector3( deltaR[DP], deltaR[dT], deltaR[dC1] );
crossProductVector3( deltaR[dU], deltaR[QD], deltaR[dC2] );
crossProductVector3(deltaR[DP], deltaR[dT], deltaR[dC1] );
crossProductVector3(deltaR[dU], deltaR[QD], deltaR[dC2] );
crossProductVector3( deltaR[dT], deltaR[CP], deltaR[dD1] );
crossProductVector3( deltaR[QC], deltaR[dU], deltaR[dD2] );
crossProductVector3(deltaR[dT], deltaR[CP], deltaR[dD1] );
crossProductVector3(deltaR[QC], deltaR[dU], deltaR[dD2] );
crossProductVector3( deltaR[BD], deltaR[dP], d[A] );
crossProductVector3( deltaR[dP], deltaR[AD], d[B] );
crossProductVector3(deltaR[BD], deltaR[dP], d[A] );
crossProductVector3(deltaR[dP], deltaR[AD], d[B] );
crossProductVector3( deltaR[FC], deltaR[dQ], d[E] );
crossProductVector3( deltaR[dQ], deltaR[EC], d[F] );
crossProductVector3(deltaR[FC], deltaR[dQ], d[E] );
crossProductVector3(deltaR[dQ], deltaR[EC], d[F] );
for( int ii = 0; ii < 3; ii++ ){
for (int ii = 0; ii < 3; ii++) {
d[C][ii] = deltaR[dC1][ii] + deltaR[dC2][ii] + deltaR[dP][ii] - d[E][ii] - d[F][ii];
d[D][ii] = deltaR[dD1][ii] + deltaR[dD2][ii] + deltaR[dQ][ii] - d[A][ii] - d[B][ii];
}
......@@ -211,69 +211,69 @@ static void computeAmoebaPiTorsionForce(int bondIndex, std::vector<Vec3>& posit
return;
}
static void computeAmoebaPiTorsionForces( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) {
static void computeAmoebaPiTorsionForces(Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log) {
// get positions and zero forces
State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() );
expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){
for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
}
// calculates forces/energy
*expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaPiTorsionForce.getNumPiTorsions(); ii++ ){
computeAmoebaPiTorsionForce(ii, positions, amoebaPiTorsionForce, expectedForces, expectedEnergy, log );
for (int ii = 0; ii < amoebaPiTorsionForce.getNumPiTorsions(); ii++) {
computeAmoebaPiTorsionForce(ii, positions, amoebaPiTorsionForce, expectedForces, expectedEnergy, log);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaPiTorsionForces: expected energy=%14.7e\n", *expectedEnergy);
for (unsigned int ii = 0; ii < positions.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
return;
}
void compareWithExpectedForceAndEnergy( Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
double tolerance, const std::string& idString, FILE* log) {
void compareWithExpectedForceAndEnergy(Context& context, AmoebaPiTorsionForce& amoebaPiTorsionForce,
double tolerance, const std::string& idString, FILE* log) {
std::vector<Vec3> expectedForces;
double expectedEnergy;
computeAmoebaPiTorsionForces( context, amoebaPiTorsionForce, expectedForces, &expectedEnergy, log );
computeAmoebaPiTorsionForces(context, amoebaPiTorsionForce, expectedForces, &expectedEnergy, log);
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaPiTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaPiTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy());
for (unsigned int ii = 0; ii < forces.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
}
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
}
void testOnePiTorsion( FILE* log ) {
void testOnePiTorsion(FILE* log) {
System system;
int numberOfParticles = 6;
for( int ii = 0; ii < numberOfParticles; ii++ ){
for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0);
}
......@@ -282,25 +282,25 @@ void testOnePiTorsion( FILE* log ) {
AmoebaPiTorsionForce* amoebaPiTorsionForce = new AmoebaPiTorsionForce();
double kTorsion = 6.85;
amoebaPiTorsionForce->addPiTorsion(0, 1, 2, 3, 4, 5, kTorsion );
amoebaPiTorsionForce->addPiTorsion(0, 1, 2, 3, 4, 5, kTorsion);
system.addForce(amoebaPiTorsionForce);
ASSERT(!amoebaPiTorsionForce->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, Platform::getPlatformByName( "Reference"));
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 0.262660000E+02, 0.254130000E+02, 0.284200000E+01 );
positions[1] = Vec3( 0.278860000E+02, 0.264630000E+02, 0.426300000E+01 );
positions[2] = Vec3( 0.269130000E+02, 0.266390000E+02, 0.353100000E+01 );
positions[0] = Vec3(0.262660000E+02, 0.254130000E+02, 0.284200000E+01);
positions[1] = Vec3(0.278860000E+02, 0.264630000E+02, 0.426300000E+01);
positions[2] = Vec3(0.269130000E+02, 0.266390000E+02, 0.353100000E+01);
positions[3] = Vec3( 0.245568230E+02, 0.250215290E+02, 0.796852800E+01 );
positions[4] = Vec3( 0.261000000E+02, 0.292530000E+02, 0.520200000E+01 );
positions[5] = Vec3( 0.254124630E+02, 0.234691880E+02, 0.773335400E+01 );
positions[3] = Vec3(0.245568230E+02, 0.250215290E+02, 0.796852800E+01);
positions[4] = Vec3(0.261000000E+02, 0.292530000E+02, 0.520200000E+01);
positions[5] = Vec3(0.254124630E+02, 0.234691880E+02, 0.773335400E+01);
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log );
compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log);
// Try changing the torsion parameters and make sure it's still correct.
......@@ -308,29 +308,29 @@ void testOnePiTorsion( FILE* log ) {
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log );
compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log);
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaPiTorsionForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log );
compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion", log);
}
int main( int numberOfArguments, char* argv[] ) {
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestReferenceAmoebaPiTorsionForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories();
FILE* log = NULL;
//FILE* log = stderr;
//FILE* log = fopen( "AmoebaPiTorsionForce1.log", "w" );;
//FILE* log = fopen("AmoebaPiTorsionForce1.log", "w");;
testOnePiTorsion( log );
testOnePiTorsion(log);
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
if (log && log != stderr)
(void) fclose(log);
#endif
}
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaStretchBendForce.cpp
View file @ 41cd79a5
......@@ -65,7 +65,7 @@ const double TOL = 1e-4;
--------------------------------------------------------------------------------------- */
static void crossProductVector3( double* vectorX, double* vectorY, double* vectorZ ){
static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
vectorZ[0] = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
vectorZ[1] = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
......@@ -74,13 +74,13 @@ static void crossProductVector3( double* vectorX, double* vectorY, double* vecto
return;
}
static double dotVector3( double* vectorX, double* vectorY ){
static double dotVector3(double* vectorX, double* vectorY) {
return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
}
static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& positions, AmoebaStretchBendForce& amoebaStretchBendForce,
std::vector<Vec3>& forces, double* energy, FILE* log ) {
std::vector<Vec3>& forces, double* energy, FILE* log) {
int particle1, particle2, particle3;
double abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend;
......@@ -88,10 +88,10 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend);
angleStretchBend *= RADIAN;
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForce: bond %d [%d %d %d] ab=%10.3e cb=%10.3e angle=%10.3e k1=%10.3e k2=%10.3e\n",
bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend );
(void) fflush( log );
if (log) {
(void) fprintf(log, "computeAmoebaStretchBendForce: bond %d [%d %d %d] ab=%10.3e cb=%10.3e angle=%10.3e k1=%10.3e k2=%10.3e\n",
bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend);
(void) fflush(log);
}
#endif
......@@ -106,31 +106,31 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
double deltaR[LastDeltaIndex][3];
double rAB2 = 0.0;
double rCB2 = 0.0;
for( int ii = 0; ii < 3; ii++ ){
for (int ii = 0; ii < 3; ii++) {
deltaR[AB][ii] = positions[particle1][ii] - positions[particle2][ii];
rAB2 += deltaR[AB][ii]*deltaR[AB][ii];
deltaR[CB][ii] = positions[particle3][ii] - positions[particle2][ii];
rCB2 += deltaR[CB][ii]*deltaR[CB][ii];
}
double rAB = sqrt( rAB2 );
double rCB = sqrt( rCB2 );
double rAB = sqrt(rAB2);
double rCB = sqrt(rCB2);
crossProductVector3( deltaR[CB], deltaR[AB], deltaR[CBxAB] );
double rP = dotVector3( deltaR[CBxAB], deltaR[CBxAB] );
rP = sqrt( rP );
crossProductVector3(deltaR[CB], deltaR[AB], deltaR[CBxAB]);
double rP = dotVector3(deltaR[CBxAB], deltaR[CBxAB]);
rP = sqrt(rP);
if( rP <= 0.0 ){
if (rP <= 0.0) {
return;
}
double dot = dotVector3( deltaR[CB], deltaR[AB] );
double dot = dotVector3(deltaR[CB], deltaR[AB]);
double cosine = dot/(rAB*rCB);
double angle;
if( cosine >= 1.0 ){
if (cosine >= 1.0) {
angle = 0.0;
}
else if( cosine <= -1.0 ){
else if (cosine <= -1.0) {
angle = PI_M;
}
else {
......@@ -142,9 +142,9 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
// P = CBxAB
crossProductVector3( deltaR[AB], deltaR[CBxAB], deltaR[ABxP] );
crossProductVector3( deltaR[CB], deltaR[CBxAB], deltaR[CBxP] );
for( int ii = 0; ii < 3; ii++ ){
crossProductVector3(deltaR[AB], deltaR[CBxAB], deltaR[ABxP]);
crossProductVector3(deltaR[CB], deltaR[CBxAB], deltaR[CBxP]);
for (int ii = 0; ii < 3; ii++) {
deltaR[ABxP][ii] *= termA;
deltaR[CBxP][ii] *= termC;
}
......@@ -162,14 +162,14 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
// forces
// calculate forces for atoms a, b, c
// the force for b is then -( a + c)
// the force for b is then -(a + c)
double subForce[LastAtomIndex][3];
double dt = angle - angleStretchBend;
for( int jj = 0; jj < 3; jj++ ){
for (int jj = 0; jj < 3; jj++) {
subForce[A][jj] = kStretchBend*dt*termA*deltaR[AB][jj] + drkk*deltaR[ABxP][jj];
subForce[C][jj] = k2StretchBend*dt*termC*deltaR[CB][jj] + drkk*deltaR[CBxP][jj];
subForce[B][jj] = -( subForce[A][jj] + subForce[C][jj] );
subForce[B][jj] = -(subForce[A][jj] + subForce[C][jj]);
}
// ---------------------------------------------------------------------------------------
......@@ -190,78 +190,78 @@ static void computeAmoebaStretchBendForce(int bondIndex, std::vector<Vec3>& pos
*energy += dt*drkk;
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForce: angle=%10.3e dt=%10.3e dr=%10.3e\n", angle, dt, dr );
(void) fflush( log );
if (log) {
(void) fprintf(log, "computeAmoebaStretchBendForce: angle=%10.3e dt=%10.3e dr=%10.3e\n", angle, dt, dr);
(void) fflush(log);
}
#endif
return;
}
static void computeAmoebaStretchBendForces( Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log ) {
static void computeAmoebaStretchBendForces(Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log) {
// get positions and zero forces
State state = context.getState(State::Positions);
std::vector<Vec3> positions = state.getPositions();
expectedForces.resize( positions.size() );
expectedForces.resize(positions.size());
for( unsigned int ii = 0; ii < expectedForces.size(); ii++ ){
for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
}
// calculates forces/energy
*expectedEnergy = 0.0;
for( int ii = 0; ii < amoebaStretchBendForce.getNumStretchBends(); ii++ ){
computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy, log );
for (int ii = 0; ii < amoebaStretchBendForce.getNumStretchBends(); ii++) {
computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy, log);
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForces: expected energy=%14.7e\n", *expectedEnergy );
for( unsigned int ii = 0; ii < positions.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaStretchBendForces: expected energy=%14.7e\n", *expectedEnergy);
for (unsigned int ii = 0; ii < positions.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
return;
}
void compareWithExpectedForceAndEnergy( Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
double tolerance, const std::string& idString, FILE* log) {
void compareWithExpectedForceAndEnergy(Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
double tolerance, const std::string& idString, FILE* log) {
std::vector<Vec3> expectedForces;
double expectedEnergy;
computeAmoebaStretchBendForces( context, amoebaStretchBendForce, expectedForces, &expectedEnergy, log );
computeAmoebaStretchBendForces(context, amoebaStretchBendForce, expectedForces, &expectedEnergy, log);
State state = context.getState(State::Forces | State::Energy);
const std::vector<Vec3> forces = state.getForces();
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaStretchBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
if (log) {
(void) fprintf(log, "computeAmoebaStretchBendForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy());
for (unsigned int ii = 0; ii < forces.size(); ii++) {
(void) fprintf(log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2]);
}
(void) fflush( log );
(void) fflush(log);
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
}
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
}
void testOneStretchBend( FILE* log ) {
void testOneStretchBend(FILE* log) {
System system;
int numberOfParticles = 3;
for( int ii = 0; ii < numberOfParticles; ii++ ){
for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0);
}
......@@ -275,21 +275,21 @@ void testOneStretchBend( FILE* log ) {
//double kStretchBend = 0.750491578E-01;
double kStretchBend = 1.0;
amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend );
amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend);
system.addForce(amoebaStretchBendForce);
ASSERT(!amoebaStretchBendForce->usesPeriodicBoundaryConditions());
ASSERT(!system.usesPeriodicBoundaryConditions());
Context context(system, integrator, Platform::getPlatformByName( "Reference"));
Context context(system, integrator, Platform::getPlatformByName("Reference"));
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 0.262660000E+02, 0.254130000E+02, 0.284200000E+01 );
positions[1] = Vec3( 0.273400000E+02, 0.244300000E+02, 0.261400000E+01 );
positions[2] = Vec3( 0.269573220E+02, 0.236108860E+02, 0.216376800E+01 );
positions[0] = Vec3(0.262660000E+02, 0.254130000E+02, 0.284200000E+01);
positions[1] = Vec3(0.273400000E+02, 0.244300000E+02, 0.261400000E+01);
positions[2] = Vec3(0.269573220E+02, 0.236108860E+02, 0.216376800E+01);
context.setPositions(positions);
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log );
compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log);
// Try changing the stretch-bend parameters and make sure it's still correct.
......@@ -297,17 +297,17 @@ void testOneStretchBend( FILE* log ) {
bool exceptionThrown = false;
try {
// This should throw an exception.
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log );
compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log);
}
catch (std::exception ex) {
exceptionThrown = true;
}
ASSERT(exceptionThrown);
amoebaStretchBendForce->updateParametersInContext(context);
compareWithExpectedForceAndEnergy( context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log );
compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend", log);
}
int main( int numberOfArguments, char* argv[] ) {
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestReferenceAmoebaStretchBendForce running test..." << std::endl;
......@@ -315,11 +315,11 @@ int main( int numberOfArguments, char* argv[] ) {
FILE* log = NULL;
//FILE* log = stderr;
//FILE* log = fopen( "AmoebaStretchBendForce1.log", "w" );;
testOneStretchBend( log );
//FILE* log = fopen("AmoebaStretchBendForce1.log", "w");;
testOneStretchBend(log);
#ifdef AMOEBA_DEBUG
if( log && log != stderr )
(void) fclose( log );
if (log && log != stderr)
(void) fclose(log);
#endif
}
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaTorsionTorsionForce.cpp
View file @ 41cd79a5
......@@ -2559,20 +2559,20 @@ static double grid[4][625][6] = {
int elementCount = (includeDerivs ? 6 : 3);
std::vector<TorsionTorsionGrid> grids(4);
for( int ii = 0; ii < 4; ii++ ){
grids[ii].resize( 25 );
for( int jj = 0; jj < 25; jj++ ){
for (int ii = 0; ii < 4; ii++) {
grids[ii].resize(25);
for (int jj = 0; jj < 25; jj++) {
grids[ii][jj].resize(25);
for( int kk = 0; kk < 25; kk++ ){
for (int kk = 0; kk < 25; kk++) {
grids[ii][jj][kk].resize(elementCount);
}
}
int index = 0;
for( int jj = 0; jj < 25; jj++ ){
for( int kk = 0; kk < 25; kk++ ){
for (int jj = 0; jj < 25; jj++) {
for (int kk = 0; kk < 25; kk++) {
int jjIndex = static_cast<int>(((grid[ii][index][0] + 180.0)/15.0)+1.0e-05);
int kkIndex = static_cast<int>(((grid[ii][index][1] + 180.0)/15.0)+1.0e-05);
for( int ll = 0; ll < elementCount; ll++ ){
for (int ll = 0; ll < elementCount; ll++) {
grids[ii][kk][jj][ll] = grid[ii][index][ll];
}
index++;
......@@ -2586,7 +2586,7 @@ void testTorsionTorsion(int systemId, bool includeDerivs) {
System system;
int numberOfParticles = 6;
for( int ii = 0; ii < numberOfParticles; ii++ ){
for (int ii = 0; ii < numberOfParticles; ii++) {
system.addParticle(1.0);
}
LangevinIntegrator integrator(0.0, 0.1, 0.01);
......@@ -2598,48 +2598,48 @@ void testTorsionTorsion(int systemId, bool includeDerivs) {
std::vector<Vec3> positions(numberOfParticles);
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy;
if( systemId == 0 ){
if (systemId == 0) {
// villin: 2 19 20 21 38 25 grid=2
chiralCheckAtomIndex = 5;
gridIndex = 2;
positions[0] = Vec3( -0.422792800E+01, -0.110605910E+02, -0.508156700E+01 );
positions[1] = Vec3( -0.447153100E+01, -0.978627900E+01, -0.466405800E+01 );
positions[2] = Vec3( -0.531878400E+01, -0.940508600E+01, -0.352283100E+01 );
positions[3] = Vec3( -0.679606000E+01, -0.974353100E+01, -0.382975700E+01 );
positions[4] = Vec3( -0.760612300E+01, -0.992590200E+01, -0.275088400E+01 );
positions[5] = Vec3( -0.516893900E+01, -0.788347000E+01, -0.316943000E+01 );
positions[0] = Vec3(-0.422792800E+01, -0.110605910E+02, -0.508156700E+01);
positions[1] = Vec3(-0.447153100E+01, -0.978627900E+01, -0.466405800E+01);
positions[2] = Vec3(-0.531878400E+01, -0.940508600E+01, -0.352283100E+01);
positions[3] = Vec3(-0.679606000E+01, -0.974353100E+01, -0.382975700E+01);
positions[4] = Vec3(-0.760612300E+01, -0.992590200E+01, -0.275088400E+01);
positions[5] = Vec3(-0.516893900E+01, -0.788347000E+01, -0.316943000E+01);
expectedForces[0] = Vec3( 0.906091624E+00, -0.529814945E-01, 0.690384140E+00 );
expectedForces[1] = Vec3( -0.124550232E+01, -0.999341692E+00, -0.590867130E+00 );
expectedForces[2] = Vec3( 0.534419689E+00, 0.612404926E-01, 0.547380310E-01 );
expectedForces[3] = Vec3( -5.732010432E-01, 2.645718463E+00, -1.585204274E-01 );
expectedForces[4] = Vec3( 3.781920539E-01, -1.654635768E+00, 4.265386268E-03 );
expectedForces[5] = Vec3( 0.0, 0.0, 0.0 );
expectedForces[0] = Vec3( 0.906091624E+00, -0.529814945E-01, 0.690384140E+00);
expectedForces[1] = Vec3(-0.124550232E+01, -0.999341692E+00, -0.590867130E+00);
expectedForces[2] = Vec3( 0.534419689E+00, 0.612404926E-01, 0.547380310E-01);
expectedForces[3] = Vec3(-5.732010432E-01, 2.645718463E+00, -1.585204274E-01);
expectedForces[4] = Vec3( 3.781920539E-01, -1.654635768E+00, 4.265386268E-03);
expectedForces[5] = Vec3( 0.0, 0.0, 0.0);
expectedEnergy = -2.699654759E+00;
}
else if( systemId == 1 ){
else if (systemId == 1) {
// villin: 158 176 177 178 183 -1 0
chiralCheckAtomIndex = -1;
gridIndex = 0;
positions[0] = Vec3( -0.105946640E+02, -0.917797000E+00, 0.105486310E+02 );
positions[1] = Vec3( -0.115059090E+02, -0.141876700E+01, 0.966933200E+01 );
positions[2] = Vec3( -0.128314660E+02, -0.876338000E+00, 0.942959800E+01 );
positions[3] = Vec3( -0.130879850E+02, -0.760280000E-01, 0.814732200E+01 );
positions[4] = Vec3( -0.120888080E+02, 0.112050000E-01, 0.722704500E+01 );
positions[5] = Vec3( 0.0, 0.0, 0.0 );
positions[0] = Vec3(-0.105946640E+02, -0.917797000E+00, 0.105486310E+02);
positions[1] = Vec3(-0.115059090E+02, -0.141876700E+01, 0.966933200E+01);
positions[2] = Vec3(-0.128314660E+02, -0.876338000E+00, 0.942959800E+01);
positions[3] = Vec3(-0.130879850E+02, -0.760280000E-01, 0.814732200E+01);
positions[4] = Vec3(-0.120888080E+02, 0.112050000E-01, 0.722704500E+01);
positions[5] = Vec3( 0.0, 0.0, 0.0);
expectedForces[0] = Vec3( 4.165851130E-01, 6.608242922E-01, -8.082168261E-01 );
expectedForces[1] = Vec3( -6.024659721E-01, -8.878744406E-01, 1.322274444E+00 );
expectedForces[2] = Vec3( 3.196925118E-02, -3.137497848E-01, -8.207984001E-01 );
expectedForces[3] = Vec3( 3.842205941E-02, 2.602732089E-01, 1.547586195E-01 );
expectedForces[4] = Vec3( 1.154895485E-01, 2.805267242E-01, 1.519821623E-01 );
expectedForces[5] = Vec3( 0.0, 0.0, 0.0 );
expectedForces[0] = Vec3( 4.165851130E-01, 6.608242922E-01, -8.082168261E-01);
expectedForces[1] = Vec3(-6.024659721E-01, -8.878744406E-01, 1.322274444E+00);
expectedForces[2] = Vec3( 3.196925118E-02, -3.137497848E-01, -8.207984001E-01);
expectedForces[3] = Vec3( 3.842205941E-02, 2.602732089E-01, 1.547586195E-01);
expectedForces[4] = Vec3( 1.154895485E-01, 2.805267242E-01, 1.519821623E-01);
expectedForces[5] = Vec3( 0.0, 0.0, 0.0);
expectedEnergy = -3.372536909E+00;
}
......@@ -2656,21 +2656,21 @@ void testTorsionTorsion(int systemId, bool includeDerivs) {
std::vector<Vec3> forces = state.getForces();
const double conversion = -1.0;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
for (unsigned int ii = 0; ii < forces.size(); ii++) {
forces[ii][0] *= conversion;
forces[ii][1] *= conversion;
forces[ii][2] *= conversion;
}
double tolerance = 1.0e-03;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
}
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
}
int main( int numberOfArguments, char* argv[] ) {
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestReferenceAmoebaTorsionTorsionForce running test..." << std::endl;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaVdwForce.cpp
View file @ 41cd79a5
This source diff could not be displayed because it is too large. You can view the blob instead.
plugins/amoeba/platforms/reference/tests/TestReferenceWcaDispersionForce.cpp
View file @ 41cd79a5
......@@ -55,21 +55,21 @@ extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();
const double TOL = 1e-4;
void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy,
const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log ) {
void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy,
const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance, FILE* log) {
// beginning of WcaDispersion setup
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC_MOD( expectedForces[ii], forces[ii], tolerance, testName );
for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC_MOD(expectedForces[ii], forces[ii], tolerance, testName);
}
ASSERT_EQUAL_TOL_MOD( expectedEnergy, energy, tolerance, testName );
ASSERT_EQUAL_TOL_MOD(expectedEnergy, energy, tolerance, testName);
}
// test Wca dispersion
void testWcaDispersionAmmonia( FILE* log ) {
void testWcaDispersionAmmonia(FILE* log) {
std::string testName = "testWcaDispersionAmmonia";
int numberOfParticles = 8;
......@@ -79,41 +79,41 @@ void testWcaDispersionAmmonia( FILE* log ) {
System system;
AmoebaWcaDispersionForce* amoebaWcaDispersionForce = new AmoebaWcaDispersionForce();;
amoebaWcaDispersionForce->setEpso( 4.6024000e-01 );
amoebaWcaDispersionForce->setEpsh( 5.6484000e-02 );
amoebaWcaDispersionForce->setRmino( 1.7025000e-01 );
amoebaWcaDispersionForce->setRminh( 1.3275000e-01 );
amoebaWcaDispersionForce->setDispoff( 2.6000000e-02 );
amoebaWcaDispersionForce->setAwater( 3.3428000e+01 );
amoebaWcaDispersionForce->setSlevy( 1.0000000e+00 );
amoebaWcaDispersionForce->setShctd( 8.1000000e-01 );
amoebaWcaDispersionForce->setEpso( 4.6024000e-01);
amoebaWcaDispersionForce->setEpsh( 5.6484000e-02);
amoebaWcaDispersionForce->setRmino( 1.7025000e-01);
amoebaWcaDispersionForce->setRminh( 1.3275000e-01);
amoebaWcaDispersionForce->setDispoff(2.6000000e-02);
amoebaWcaDispersionForce->setAwater( 3.3428000e+01);
amoebaWcaDispersionForce->setSlevy( 1.0000000e+00);
amoebaWcaDispersionForce->setShctd( 8.1000000e-01);
// addParticle: radius, epsilon
for( unsigned int ii = 0; ii < 2; ii++ ){
system.addParticle( 1.4007000e+01 );
amoebaWcaDispersionForce->addParticle( 1.8550000e-01, 4.3932000e-01 );
for (unsigned int ii = 0; ii < 2; ii++) {
system.addParticle( 1.4007000e+01);
amoebaWcaDispersionForce->addParticle( 1.8550000e-01, 4.3932000e-01);
system.addParticle( 1.0080000e+00 );
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02 );
system.addParticle( 1.0080000e+00);
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02);
system.addParticle( 1.0080000e+00 );
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02 );
system.addParticle( 1.0080000e+00);
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02);
system.addParticle( 1.0080000e+00 );
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02 );
system.addParticle( 1.0080000e+00);
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02);
}
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03 );
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02 );
positions[2] = Vec3( 2.0843610e-01, 8.0953200e-02, 3.7462200e-02 );
positions[3] = Vec3( 1.7280780e-01, 2.0730000e-04, -9.8741700e-02 );
positions[4] = Vec3( -1.6743680e-01, 1.5900000e-05, -6.6149000e-03 );
positions[5] = Vec3( -2.0428260e-01, 8.1071500e-02, 4.1343900e-02 );
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02 );
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02 );
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03);
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02);
positions[2] = Vec3( 2.0843610e-01, 8.0953200e-02, 3.7462200e-02);
positions[3] = Vec3( 1.7280780e-01, 2.0730000e-04, -9.8741700e-02);
positions[4] = Vec3( -1.6743680e-01, 1.5900000e-05, -6.6149000e-03);
positions[5] = Vec3( -2.0428260e-01, 8.1071500e-02, 4.1343900e-02);
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02);
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02);
system.addForce(amoebaWcaDispersionForce);
ASSERT(!amoebaWcaDispersionForce->usesPeriodicBoundaryConditions());
......@@ -122,7 +122,7 @@ void testWcaDispersionAmmonia( FILE* log ) {
std::string platformName;
platformName = "Reference";
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName( platformName ) );
Context context(system, integrator, Platform::getPlatformByName(platformName));
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
......@@ -134,17 +134,17 @@ void testWcaDispersionAmmonia( FILE* log ) {
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = -2.6981209e+01;
expectedForces[0] = Vec3( 4.7839388e+00, -7.3510133e-04, -5.0382764e-01 );
expectedForces[1] = Vec3( 1.4657758e+00, 1.2431003e+00, -6.7075886e-01 );
expectedForces[2] = Vec3( 1.4563936e+00, -1.2399917e+00, -6.7443841e-01 );
expectedForces[3] = Vec3( 2.1116744e+00, -2.7407512e-03, 1.3271245e+00 );
expectedForces[4] = Vec3( -4.7528440e+00, -1.5148066e-03, 1.2653813e+00 );
expectedForces[5] = Vec3( -1.1875619e+00, -1.2866678e+00, -3.9109060e-01 );
expectedForces[6] = Vec3( -2.6885679e+00, -4.3038639e-04, 3.3763583e-02 );
expectedForces[7] = Vec3( -1.1888087e+00, 1.2889802e+00, -3.8615387e-01 );
expectedForces[0] = Vec3( 4.7839388e+00, -7.3510133e-04, -5.0382764e-01);
expectedForces[1] = Vec3( 1.4657758e+00, 1.2431003e+00, -6.7075886e-01);
expectedForces[2] = Vec3( 1.4563936e+00, -1.2399917e+00, -6.7443841e-01);
expectedForces[3] = Vec3( 2.1116744e+00, -2.7407512e-03, 1.3271245e+00);
expectedForces[4] = Vec3( -4.7528440e+00, -1.5148066e-03, 1.2653813e+00);
expectedForces[5] = Vec3( -1.1875619e+00, -1.2866678e+00, -3.9109060e-01);
expectedForces[6] = Vec3( -2.6885679e+00, -4.3038639e-04, 3.3763583e-02);
expectedForces[7] = Vec3( -1.1888087e+00, 1.2889802e+00, -3.8615387e-01);
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance, log);
// Try changing the particle parameters and make sure it's still correct.
......@@ -172,7 +172,7 @@ void testWcaDispersionAmmonia( FILE* log ) {
compareForcesEnergy(testName, state1.getPotentialEnergy(), state2.getPotentialEnergy(), state1.getForces(), state2.getForces(), tolerance, log);
}
int main( int numberOfArguments, char* argv[] ) {
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestCudaAmoebaWcaDispersionForce running test..." << std::endl;
......@@ -182,7 +182,7 @@ int main( int numberOfArguments, char* argv[] ) {
// test Wca dispersion force using two ammonia molecules
testWcaDispersionAmmonia( log );
testWcaDispersionAmmonia(log);
}
......
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