Skip to content

GitLab

Commit a783b996 authored by peastman's avatar peastman
Browse files

Eliminated RealOpenMM type

parent 9500f3af
Hide whitespace changes
Inline Side-by-side
Showing with 1056 additions and 1326 deletions
platforms/reference/src/SimTKReference/ReferenceNeighborList.cpp
View file @ a783b996
......@@ -13,8 +13,8 @@ namespace OpenMM {
typedef std::vector<AtomIndex> AtomList;
// squared distance between two points
static double compPairDistanceSquared(const RealVec& pos1, const RealVec& pos2, const RealVec* periodicBoxVectors, bool usePeriodic) {
RealVec diff = pos2-pos1;
static double compPairDistanceSquared(const Vec3& pos1, const Vec3& pos2, const Vec3* periodicBoxVectors, bool usePeriodic) {
Vec3 diff = pos2-pos1;
if (usePeriodic) {
diff -= periodicBoxVectors[2]*floor(diff[2]/periodicBoxVectors[2][2]+0.5);
diff -= periodicBoxVectors[1]*floor(diff[1]/periodicBoxVectors[1][1]+0.5);
......@@ -30,7 +30,7 @@ void OPENMM_EXPORT computeNeighborListNaive(
int nAtoms,
const AtomLocationList& atomLocations,
const vector<set<int> >& exclusions,
const RealVec* periodicBoxVectors,
const Vec3* periodicBoxVectors,
bool usePeriodic,
double maxDistance,
double minDistance,
......@@ -78,13 +78,13 @@ public:
};
typedef std::pair<const RealVec*, AtomIndex> VoxelItem;
typedef std::pair<const Vec3*, AtomIndex> VoxelItem;
typedef std::vector< VoxelItem > Voxel;
class VoxelHash
{
public:
VoxelHash(double vsx, double vsy, double vsz, const RealVec* periodicBoxVectors, bool usePeriodic) :
VoxelHash(double vsx, double vsy, double vsz, const Vec3* periodicBoxVectors, bool usePeriodic) :
voxelSizeX(vsx), voxelSizeY(vsy), voxelSizeZ(vsz), periodicBoxVectors(periodicBoxVectors), usePeriodic(usePeriodic) {
if (usePeriodic) {
nx = (int) floor(periodicBoxVectors[0][0]/voxelSizeX+0.5);
......@@ -96,7 +96,7 @@ public:
}
}
void insert(const AtomIndex& item, const RealVec& location)
void insert(const AtomIndex& item, const Vec3& location)
{
VoxelIndex voxelIndex = getVoxelIndex(location);
if (voxelMap.find(voxelIndex) == voxelMap.end())
......@@ -106,8 +106,8 @@ public:
}
VoxelIndex getVoxelIndex(const RealVec& location) const {
RealVec r = location;
VoxelIndex getVoxelIndex(const Vec3& location) const {
Vec3 r = location;
if (usePeriodic) {
r -= periodicBoxVectors[2]*floor(r[2]/periodicBoxVectors[2][2]);
r -= periodicBoxVectors[1]*floor(r[1]/periodicBoxVectors[1][1]);
......@@ -138,7 +138,7 @@ public:
assert(voxelSizeZ > 0);
const AtomIndex atomI = referencePoint.second;
const RealVec& locationI = *referencePoint.first;
const Vec3& locationI = *referencePoint.first;
double maxDistanceSquared = maxDistance * maxDistance;
double minDistanceSquared = minDistance * minDistance;
......@@ -187,7 +187,7 @@ public:
for (Voxel::const_iterator itemIter = voxel.begin(); itemIter != voxel.end(); ++itemIter)
{
const AtomIndex atomJ = itemIter->second;
const RealVec& locationJ = *itemIter->first;
const Vec3& locationJ = *itemIter->first;
// Ignore self hits
if (atomI == atomJ) continue;
......@@ -211,7 +211,7 @@ public:
private:
double voxelSizeX, voxelSizeY, voxelSizeZ;
int nx, ny, nz;
const RealVec* periodicBoxVectors;
const Vec3* periodicBoxVectors;
const bool usePeriodic;
std::map<VoxelIndex, Voxel> voxelMap;
};
......@@ -223,7 +223,7 @@ void OPENMM_EXPORT computeNeighborListVoxelHash(
int nAtoms,
const AtomLocationList& atomLocations,
const vector<set<int> >& exclusions,
const RealVec* periodicBoxVectors,
const Vec3* periodicBoxVectors,
bool usePeriodic,
double maxDistance,
double minDistance,
......@@ -244,7 +244,7 @@ void OPENMM_EXPORT computeNeighborListVoxelHash(
for (AtomIndex atomJ = 0; atomJ < (AtomIndex) nAtoms; ++atomJ) // use "j", because j > i for pairs
{
// 1) Find other atoms that are close to this one
const RealVec& location = atomLocations[atomJ];
const Vec3& location = atomLocations[atomJ];
voxelHash.getNeighbors(
neighborList,
VoxelItem(&location, atomJ),
......
platforms/reference/src/SimTKReference/ReferenceObc.cpp
View file @ a783b996
......@@ -111,7 +111,7 @@ void ReferenceObc::setIncludeAceApproximation(int includeAceApproximation) {
--------------------------------------------------------------------------------------- */
vector<RealOpenMM>& ReferenceObc::getObcChain() {
vector<double>& ReferenceObc::getObcChain() {
return _obcChain;
}
......@@ -127,30 +127,19 @@ vector<RealOpenMM>& ReferenceObc::getObcChain() {
--------------------------------------------------------------------------------------- */
void ReferenceObc::computeBornRadii(const vector<RealVec>& atomCoordinates, vector<RealOpenMM>& bornRadii) {
void ReferenceObc::computeBornRadii(const vector<Vec3>& atomCoordinates, vector<double>& bornRadii) {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>(0.0);
static const RealOpenMM one = static_cast<RealOpenMM>(1.0);
static const RealOpenMM two = static_cast<RealOpenMM>(2.0);
static const RealOpenMM three = static_cast<RealOpenMM>(3.0);
static const RealOpenMM half = static_cast<RealOpenMM>(0.5);
static const RealOpenMM fourth = static_cast<RealOpenMM>(0.25);
// ---------------------------------------------------------------------------------------
ObcParameters* obcParameters = getObcParameters();
ObcParameters* obcParameters = getObcParameters();
int numberOfAtoms = obcParameters->getNumberOfAtoms();
const vector<double>& atomicRadii = obcParameters->getAtomicRadii();
const vector<double>& scaledRadiusFactor = obcParameters->getScaledRadiusFactors();
vector<double>& obcChain = getObcChain();
int numberOfAtoms = obcParameters->getNumberOfAtoms();
const vector<RealOpenMM>& atomicRadii = obcParameters->getAtomicRadii();
const vector<RealOpenMM>& scaledRadiusFactor = obcParameters->getScaledRadiusFactors();
vector<RealOpenMM>& obcChain = getObcChain();
RealOpenMM dielectricOffset = obcParameters->getDielectricOffset();
RealOpenMM alphaObc = obcParameters->getAlphaObc();
RealOpenMM betaObc = obcParameters->getBetaObc();
RealOpenMM gammaObc = obcParameters->getGammaObc();
double dielectricOffset = obcParameters->getDielectricOffset();
double alphaObc = obcParameters->getAlphaObc();
double betaObc = obcParameters->getBetaObc();
double gammaObc = obcParameters->getGammaObc();
// ---------------------------------------------------------------------------------------
......@@ -158,11 +147,11 @@ void ReferenceObc::computeBornRadii(const vector<RealVec>& atomCoordinates, vect
for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
double radiusI = atomicRadii[atomI];
double offsetRadiusI = radiusI - dielectricOffset;
RealOpenMM radiusIInverse = one/offsetRadiusI;
RealOpenMM sum = zero;
double radiusIInverse = 1.0/offsetRadiusI;
double sum = 0.0;
// HCT code
......@@ -170,38 +159,38 @@ void ReferenceObc::computeBornRadii(const vector<RealVec>& atomCoordinates, vect
if (atomJ != atomI) {
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
double deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR);
else
ReferenceForce::getDeltaR(atomCoordinates[atomI], atomCoordinates[atomJ], deltaR);
RealOpenMM r = deltaR[ReferenceForce::RIndex];
double r = deltaR[ReferenceForce::RIndex];
if (_obcParameters->getUseCutoff() && r > _obcParameters->getCutoffDistance())
continue;
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
double offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
double scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
double rScaledRadiusJ = r + scaledRadiusJ;
if (offsetRadiusI < rScaledRadiusJ) {
RealOpenMM rInverse = one/r;
RealOpenMM l_ij = offsetRadiusI > FABS(r - scaledRadiusJ) ? offsetRadiusI : FABS(r - scaledRadiusJ);
l_ij = one/l_ij;
double rInverse = 1.0/r;
double l_ij = offsetRadiusI > fabs(r - scaledRadiusJ) ? offsetRadiusI : fabs(r - scaledRadiusJ);
l_ij = 1.0/l_ij;
RealOpenMM u_ij = one/rScaledRadiusJ;
double u_ij = 1.0/rScaledRadiusJ;
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
double l_ij2 = l_ij*l_ij;
double u_ij2 = u_ij*u_ij;
RealOpenMM ratio = LN((u_ij/l_ij));
RealOpenMM term = l_ij - u_ij + fourth*r*(u_ij2 - l_ij2) + (half*rInverse*ratio) + (fourth*scaledRadiusJ*scaledRadiusJ*rInverse)*(l_ij2 - u_ij2);
double ratio = log((u_ij/l_ij));
double term = l_ij - u_ij + 0.25*r*(u_ij2 - l_ij2) + (0.5*rInverse*ratio) + (0.25*scaledRadiusJ*scaledRadiusJ*rInverse)*(l_ij2 - u_ij2);
// this case (atom i completely inside atom j) is not considered in the original paper
// Jay Ponder and the authors of Tinker recognized this and
// worked out the details
if (offsetRadiusI < (scaledRadiusJ - r)) {
term += two*(radiusIInverse - l_ij);
term += 2.0*(radiusIInverse - l_ij);
}
sum += term;
......@@ -211,15 +200,15 @@ void ReferenceObc::computeBornRadii(const vector<RealVec>& atomCoordinates, vect
// OBC-specific code (Eqs. 6-8 in paper)
sum *= half*offsetRadiusI;
RealOpenMM sum2 = sum*sum;
RealOpenMM sum3 = sum*sum2;
RealOpenMM tanhSum = TANH(alphaObc*sum - betaObc*sum2 + gammaObc*sum3);
sum *= 0.5*offsetRadiusI;
double sum2 = sum*sum;
double sum3 = sum*sum2;
double tanhSum = tanh(alphaObc*sum - betaObc*sum2 + gammaObc*sum3);
bornRadii[atomI] = one/(one/offsetRadiusI - tanhSum/radiusI);
bornRadii[atomI] = 1.0/(1.0/offsetRadiusI - tanhSum/radiusI);
obcChain[atomI] = offsetRadiusI*(alphaObc - two*betaObc*sum + three*gammaObc*sum2);
obcChain[atomI] = (one - tanhSum*tanhSum)*obcChain[atomI]/radiusI;
obcChain[atomI] = offsetRadiusI*(alphaObc - 2.0*betaObc*sum + 3.0*gammaObc*sum2);
obcChain[atomI] = (1.0 - tanhSum*tanhSum)*obcChain[atomI]/radiusI;
}
}
......@@ -236,24 +225,16 @@ void ReferenceObc::computeBornRadii(const vector<RealVec>& atomCoordinates, vect
--------------------------------------------------------------------------------------- */
void ReferenceObc::computeAceNonPolarForce(const ObcParameters* obcParameters,
const vector<RealOpenMM>& bornRadii,
RealOpenMM* energy,
vector<RealOpenMM>& forces) const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>(0.0);
static const RealOpenMM minusSix = -6.0;
static const RealOpenMM six = static_cast<RealOpenMM>(6.0);
// ---------------------------------------------------------------------------------------
const vector<double>& bornRadii,
double* energy,
vector<double>& forces) const {
// compute the nonpolar solvation via ACE approximation
const RealOpenMM probeRadius = obcParameters->getProbeRadius();
const RealOpenMM surfaceAreaFactor = obcParameters->getPi4Asolv();
const double probeRadius = obcParameters->getProbeRadius();
const double surfaceAreaFactor = obcParameters->getPi4Asolv();
const vector<RealOpenMM>& atomicRadii = obcParameters->getAtomicRadii();
const vector<double>& atomicRadii = obcParameters->getAtomicRadii();
int numberOfAtoms = obcParameters->getNumberOfAtoms();
// the original ACE equation is based on Eq.2 of
......@@ -271,12 +252,12 @@ void ReferenceObc::computeAceNonPolarForce(const ObcParameters* obcParameters,
// no paper to cite.
for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
if (bornRadii[atomI] > zero) {
RealOpenMM r = atomicRadii[atomI] + probeRadius;
RealOpenMM ratio6 = POW(atomicRadii[atomI]/bornRadii[atomI], six);
RealOpenMM saTerm = surfaceAreaFactor*r*r*ratio6;
if (bornRadii[atomI] > 0.0) {
double r = atomicRadii[atomI] + probeRadius;
double ratio6 = pow(atomicRadii[atomI]/bornRadii[atomI], 6.0);
double saTerm = surfaceAreaFactor*r*r*ratio6;
*energy += saTerm;
forces[atomI] += minusSix*saTerm/bornRadii[atomI];
forces[atomI] -= 6.0*saTerm/bornRadii[atomI];
}
}
}
......@@ -293,44 +274,33 @@ void ReferenceObc::computeAceNonPolarForce(const ObcParameters* obcParameters,
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceObc::computeBornEnergyForces(const vector<RealVec>& atomCoordinates,
const vector<RealOpenMM>& partialCharges, vector<RealVec>& inputForces) {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>(0.0);
static const RealOpenMM one = static_cast<RealOpenMM>(1.0);
static const RealOpenMM two = static_cast<RealOpenMM>(2.0);
static const RealOpenMM three = static_cast<RealOpenMM>(3.0);
static const RealOpenMM four = static_cast<RealOpenMM>(4.0);
static const RealOpenMM half = static_cast<RealOpenMM>(0.5);
static const RealOpenMM fourth = static_cast<RealOpenMM>(0.25);
static const RealOpenMM eighth = static_cast<RealOpenMM>(0.125);
double ReferenceObc::computeBornEnergyForces(const vector<Vec3>& atomCoordinates,
const vector<double>& partialCharges, vector<Vec3>& inputForces) {
// constants
const int numberOfAtoms = _obcParameters->getNumberOfAtoms();
const RealOpenMM dielectricOffset = _obcParameters->getDielectricOffset();
const RealOpenMM cutoffDistance = _obcParameters->getCutoffDistance();
const RealOpenMM soluteDielectric = _obcParameters->getSoluteDielectric();
const RealOpenMM solventDielectric = _obcParameters->getSolventDielectric();
RealOpenMM preFactor;
if (soluteDielectric != zero && solventDielectric != zero)
preFactor = two*_obcParameters->getElectricConstant()*((one/soluteDielectric) - (one/solventDielectric));
const double dielectricOffset = _obcParameters->getDielectricOffset();
const double cutoffDistance = _obcParameters->getCutoffDistance();
const double soluteDielectric = _obcParameters->getSoluteDielectric();
const double solventDielectric = _obcParameters->getSolventDielectric();
double preFactor;
if (soluteDielectric != 0.0 && solventDielectric != 0.0)
preFactor = 2.0*_obcParameters->getElectricConstant()*((1.0/soluteDielectric) - (1.0/solventDielectric));
else
preFactor = zero;
preFactor = 0.0;
// ---------------------------------------------------------------------------------------
// compute Born radii
vector<RealOpenMM> bornRadii(numberOfAtoms);
vector<double> bornRadii(numberOfAtoms);
computeBornRadii(atomCoordinates, bornRadii);
// set energy/forces to zero
RealOpenMM obcEnergy = zero;
vector<RealOpenMM> bornForces(numberOfAtoms, 0.0);
double obcEnergy = 0.0;
vector<double> bornForces(numberOfAtoms, 0.0);
// ---------------------------------------------------------------------------------------
......@@ -346,10 +316,10 @@ RealOpenMM ReferenceObc::computeBornEnergyForces(const vector<RealVec>& atomCoor
for (int atomI = 0; atomI < numberOfAtoms; atomI++) {
RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
double partialChargeI = preFactor*partialCharges[atomI];
for (int atomJ = atomI; atomJ < numberOfAtoms; atomJ++) {
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
double deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR);
else
......@@ -357,24 +327,24 @@ RealOpenMM ReferenceObc::computeBornEnergyForces(const vector<RealVec>& atomCoor
if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > cutoffDistance)
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
double r2 = deltaR[ReferenceForce::R2Index];
double deltaX = deltaR[ReferenceForce::XIndex];
double deltaY = deltaR[ReferenceForce::YIndex];
double deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
RealOpenMM D_ij = r2/(four*alpha2_ij);
double alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
double D_ij = r2/(4.0*alpha2_ij);
RealOpenMM expTerm = EXP(-D_ij);
RealOpenMM denominator2 = r2 + alpha2_ij*expTerm;
RealOpenMM denominator = SQRT(denominator2);
double expTerm = exp(-D_ij);
double denominator2 = r2 + alpha2_ij*expTerm;
double denominator = sqrt(denominator2);
RealOpenMM Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
RealOpenMM dGpol_dr = -Gpol*(one - fourth*expTerm)/denominator2;
double Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
double dGpol_dr = -Gpol*(1.0 - 0.25*expTerm)/denominator2;
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*(one + D_ij)/denominator2;
double dGpol_dalpha2_ij = -0.5*Gpol*expTerm*(1.0 + D_ij)/denominator2;
RealOpenMM energy = Gpol;
double energy = Gpol;
if (atomI != atomJ) {
......@@ -396,7 +366,7 @@ RealOpenMM ReferenceObc::computeBornEnergyForces(const vector<RealVec>& atomCoor
inputForces[atomJ][2] -= deltaZ;
} else {
energy *= half;
energy *= 0.5;
}
obcEnergy += energy;
......@@ -409,13 +379,13 @@ RealOpenMM ReferenceObc::computeBornEnergyForces(const vector<RealVec>& atomCoor
// second main loop
const vector<RealOpenMM>& obcChain = getObcChain();
const vector<RealOpenMM>& atomicRadii = _obcParameters->getAtomicRadii();
const vector<double>& obcChain = getObcChain();
const vector<double>& atomicRadii = _obcParameters->getAtomicRadii();
const RealOpenMM alphaObc = _obcParameters->getAlphaObc();
const RealOpenMM betaObc = _obcParameters->getBetaObc();
const RealOpenMM gammaObc = _obcParameters->getGammaObc();
const vector<RealOpenMM>& scaledRadiusFactor = _obcParameters->getScaledRadiusFactors();
const double alphaObc = _obcParameters->getAlphaObc();
const double betaObc = _obcParameters->getBetaObc();
const double gammaObc = _obcParameters->getGammaObc();
const vector<double>& scaledRadiusFactor = _obcParameters->getScaledRadiusFactors();
// compute factor that depends only on the outer loop index
......@@ -427,14 +397,14 @@ RealOpenMM ReferenceObc::computeBornEnergyForces(const vector<RealVec>& atomCoor
// radius w/ dielectric offset applied
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
double radiusI = atomicRadii[atomI];
double offsetRadiusI = radiusI - dielectricOffset;
for (int atomJ = 0; atomJ < numberOfAtoms; atomJ++) {
if (atomJ != atomI) {
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
double deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomI], atomCoordinates[atomJ], _obcParameters->getPeriodicBox(), deltaR);
else
......@@ -442,39 +412,39 @@ RealOpenMM ReferenceObc::computeBornEnergyForces(const vector<RealVec>& atomCoor
if (_obcParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > cutoffDistance)
continue;
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM r = deltaR[ReferenceForce::RIndex];
double deltaX = deltaR[ReferenceForce::XIndex];
double deltaY = deltaR[ReferenceForce::YIndex];
double deltaZ = deltaR[ReferenceForce::ZIndex];
double r = deltaR[ReferenceForce::RIndex];
// radius w/ dielectric offset applied
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
double offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM scaledRadiusJ2 = scaledRadiusJ*scaledRadiusJ;
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
double scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
double scaledRadiusJ2 = scaledRadiusJ*scaledRadiusJ;
double rScaledRadiusJ = r + scaledRadiusJ;
// dL/dr & dU/dr are zero (this can be shown analytically)
// removed from calculation
if (offsetRadiusI < rScaledRadiusJ) {
RealOpenMM l_ij = offsetRadiusI > FABS(r - scaledRadiusJ) ? offsetRadiusI : FABS(r - scaledRadiusJ);
l_ij = one/l_ij;
double l_ij = offsetRadiusI > fabs(r - scaledRadiusJ) ? offsetRadiusI : fabs(r - scaledRadiusJ);
l_ij = 1.0/l_ij;
RealOpenMM u_ij = one/rScaledRadiusJ;
double u_ij = 1.0/rScaledRadiusJ;
RealOpenMM l_ij2 = l_ij*l_ij;
double l_ij2 = l_ij*l_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
double u_ij2 = u_ij*u_ij;
RealOpenMM rInverse = one/r;
RealOpenMM r2Inverse = rInverse*rInverse;
double rInverse = 1.0/r;
double r2Inverse = rInverse*rInverse;
RealOpenMM t3 = eighth*(one + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + fourth*LN(u_ij/l_ij)*r2Inverse;
double t3 = 0.125*(1.0 + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + 0.25*log(u_ij/l_ij)*r2Inverse;
RealOpenMM de = bornForces[atomI]*t3*rInverse;
double de = bornForces[atomI]*t3*rInverse;
deltaX *= de;
deltaY *= de;
......
platforms/reference/src/SimTKReference/ReferencePME.cpp
View file @ a783b996
......@@ -35,6 +35,7 @@
#include "ReferencePME.h"
#include "fftpack.h"
#include "SimTKOpenMMRealType.h"
using std::vector;
......@@ -45,7 +46,7 @@ namespace OpenMM {
struct pme
{
int natoms;
RealOpenMM ewaldcoeff;
double ewaldcoeff;
t_complex * grid; /* Memory for the grid we spread charges on.
* Element (i,j,k) is accessed as:
......@@ -57,9 +58,9 @@ struct pme
int order; /* PME interpolation order. Almost always 4 */
/* Data for bspline interpolation, see the Essman PME paper */
RealOpenMM * bsplines_moduli[3]; /* 3 pointers, to x/y/z bspline moduli, each of length ngrid[x/y/z] */
RealOpenMM * bsplines_theta[3]; /* each of x/y/z has length order*natoms */
RealOpenMM * bsplines_dtheta[3]; /* each of x/y/z has length order*natoms */
double * bsplines_moduli[3]; /* 3 pointers, to x/y/z bspline moduli, each of length ngrid[x/y/z] */
double * bsplines_theta[3]; /* each of x/y/z has length order*natoms */
double * bsplines_dtheta[3]; /* each of x/y/z has length order*natoms */
ivec * particleindex; /* Array of length natoms. Each element is
* an ivec (3 ints) that specify the grid
......@@ -85,7 +86,7 @@ struct pme
* the central cell, i.e., in this case we would assume all coordinates fall in -10 nm < x,y,z < 20 nm.
*/
RealOpenMM epsilon_r; /* Dielectric coefficient to use, typically 1.0 */
double epsilon_r; /* Dielectric coefficient to use, typically 1.0 */
};
......@@ -101,25 +102,25 @@ pme_calculate_bsplines_moduli(pme_t pme)
int i,j,k,l,d;
int order;
int ndata;
RealOpenMM * data;
RealOpenMM * ddata;
RealOpenMM * bsplines_data;
RealOpenMM div;
RealOpenMM sc,ss,arg;
double * data;
double * ddata;
double * bsplines_data;
double div;
double sc,ss,arg;
nmax = 0;
for (d=0;d<3;d++)
{
nmax = (pme->ngrid[d] > nmax) ? pme->ngrid[d] : nmax;
pme->bsplines_moduli[d] = (RealOpenMM *) malloc(sizeof(RealOpenMM)*pme->ngrid[d]);
pme->bsplines_moduli[d] = (double *) malloc(sizeof(double)*pme->ngrid[d]);
}
order = pme->order;
/* temp storage in this routine */
data = (RealOpenMM *) malloc(sizeof(RealOpenMM)*order);
ddata = (RealOpenMM *) malloc(sizeof(RealOpenMM)*order);
bsplines_data = (RealOpenMM *) malloc(sizeof(RealOpenMM)*nmax);
data = (double *) malloc(sizeof(double)*order);
ddata = (double *) malloc(sizeof(double)*order);
bsplines_data = (double *) malloc(sizeof(double)*nmax);
data[order-1]=0;
data[1]=0;
......@@ -127,7 +128,7 @@ pme_calculate_bsplines_moduli(pme_t pme)
for (k=3;k<order;k++)
{
div=(RealOpenMM) (1.0/(k-1.0));
div=1.0/(k-1.0);
data[k-1]=0;
for (l=1;l<(k-1);l++)
{
......@@ -143,7 +144,7 @@ pme_calculate_bsplines_moduli(pme_t pme)
ddata[k]=data[k-1]-data[k];
}
div=(RealOpenMM) (1.0/(order-1));
div=1.0/(order-1);
data[order-1]=0;
for (l=1;l<(order-1);l++)
......@@ -170,7 +171,7 @@ pme_calculate_bsplines_moduli(pme_t pme)
sc=ss=0;
for (j=0;j<ndata;j++)
{
arg=(RealOpenMM) ((2.0*M_PI*i*j)/ndata);
arg=(2.0*M_PI*i*j)/ndata;
sc+=bsplines_data[j]*cos(arg);
ss+=bsplines_data[j]*sin(arg);
}
......@@ -192,25 +193,25 @@ pme_calculate_bsplines_moduli(pme_t pme)
}
static void invert_box_vectors(const RealVec boxVectors[3], RealVec recipBoxVectors[3])
static void invert_box_vectors(const Vec3 boxVectors[3], Vec3 recipBoxVectors[3])
{
RealOpenMM determinant = boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2];
double determinant = boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2];
assert(determinant > 0);
RealOpenMM scale = 1.0/determinant;
recipBoxVectors[0] = RealVec(boxVectors[1][1]*boxVectors[2][2], 0, 0)*scale;
recipBoxVectors[1] = RealVec(-boxVectors[1][0]*boxVectors[2][2], boxVectors[0][0]*boxVectors[2][2], 0)*scale;
recipBoxVectors[2] = RealVec(boxVectors[1][0]*boxVectors[2][1]-boxVectors[1][1]*boxVectors[2][0], -boxVectors[0][0]*boxVectors[2][1], boxVectors[0][0]*boxVectors[1][1])*scale;
double scale = 1.0/determinant;
recipBoxVectors[0] = Vec3(boxVectors[1][1]*boxVectors[2][2], 0, 0)*scale;
recipBoxVectors[1] = Vec3(-boxVectors[1][0]*boxVectors[2][2], boxVectors[0][0]*boxVectors[2][2], 0)*scale;
recipBoxVectors[2] = Vec3(boxVectors[1][0]*boxVectors[2][1]-boxVectors[1][1]*boxVectors[2][0], -boxVectors[0][0]*boxVectors[2][1], boxVectors[0][0]*boxVectors[1][1])*scale;
}
static void
pme_update_grid_index_and_fraction(pme_t pme,
const vector<RealVec>& atomCoordinates,
const RealVec periodicBoxVectors[3],
const RealVec recipBoxVectors[3])
const vector<Vec3>& atomCoordinates,
const Vec3 periodicBoxVectors[3],
const Vec3 recipBoxVectors[3])
{
int i;
int d;
RealOpenMM t;
double t;
int ti;
for (i=0;i<pme->natoms;i++)
......@@ -251,7 +252,7 @@ pme_update_grid_index_and_fraction(pme_t pme,
* numerical problems, so this shouldnt cause any problems.
* (And, by adding 100.0 box lengths, we would lose a bit of numerical accuracy here!)
*/
RealVec coord = atomCoordinates[i];
Vec3 coord = atomCoordinates[i];
for (d=0;d<3;d++)
{
t = coord[0]*recipBoxVectors[0][d]+coord[1]*recipBoxVectors[1][d]+coord[2]*recipBoxVectors[2][d];
......@@ -275,9 +276,9 @@ pme_update_bsplines(pme_t pme)
{
int i,j,k,l;
int order;
RealOpenMM dr,div;
RealOpenMM * data;
RealOpenMM * ddata;
double dr,div;
double * data;
double * ddata;
order = pme->order;
......@@ -296,7 +297,7 @@ pme_update_bsplines(pme_t pme)
for (k=3; k<order; k++)
{
div = (RealOpenMM) (1.0/(k-1.0));
div = 1.0/(k-1.0);
data[k-1] = div*dr*data[k-2];
for (l=1; l<(k-1); l++)
{
......@@ -313,7 +314,7 @@ pme_update_bsplines(pme_t pme)
ddata[k] = data[k-1]-data[k];
}
div = (RealOpenMM) (1.0/(order-1));
div = 1.0/(order-1);
data[order-1] = div*dr*data[order-2];
for (l=1; l<(order-1); l++)
......@@ -327,7 +328,7 @@ pme_update_bsplines(pme_t pme)
static void
pme_grid_spread_charge(pme_t pme, const vector<RealOpenMM>& charges)
pme_grid_spread_charge(pme_t pme, const vector<double>& charges)
{
int order;
int i;
......@@ -335,10 +336,10 @@ pme_grid_spread_charge(pme_t pme, const vector<RealOpenMM>& charges)
int x0index,y0index,z0index;
int xindex,yindex,zindex;
int index;
RealOpenMM q;
RealOpenMM * thetax;
RealOpenMM * thetay;
RealOpenMM * thetaz;
double q;
double * thetax;
double * thetay;
double * thetaz;
order = pme->order;
......@@ -407,24 +408,24 @@ pme_grid_spread_charge(pme_t pme, const vector<RealOpenMM>& charges)
static void
pme_reciprocal_convolution(pme_t pme,
const RealVec periodicBoxVectors[3],
const RealVec recipBoxVectors[3],
RealOpenMM * energy)
const Vec3 periodicBoxVectors[3],
const Vec3 recipBoxVectors[3],
double * energy)
{
int kx,ky,kz;
int nx,ny,nz;
RealOpenMM mx,my,mz;
RealOpenMM mhx,mhy,mhz,m2;
RealOpenMM one_4pi_eps;
RealOpenMM virxx,virxy,virxz,viryy,viryz,virzz;
RealOpenMM bx,by,bz;
RealOpenMM d1,d2;
RealOpenMM eterm,vfactor,struct2,ets2;
RealOpenMM esum;
RealOpenMM factor;
RealOpenMM denom;
RealOpenMM boxfactor;
RealOpenMM maxkx,maxky,maxkz;
double mx,my,mz;
double mhx,mhy,mhz,m2;
double one_4pi_eps;
double virxx,virxy,virxz,viryy,viryz,virzz;
double bx,by,bz;
double d1,d2;
double eterm,vfactor,struct2,ets2;
double esum;
double factor;
double denom;
double boxfactor;
double maxkx,maxky,maxkz;
t_complex *ptr;
......@@ -432,9 +433,9 @@ pme_reciprocal_convolution(pme_t pme,
ny = pme->ngrid[1];
nz = pme->ngrid[2];
one_4pi_eps = (RealOpenMM) (ONE_4PI_EPS0/pme->epsilon_r);
factor = (RealOpenMM) (M_PI*M_PI/(pme->ewaldcoeff*pme->ewaldcoeff));
boxfactor = (RealOpenMM) (M_PI*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2]);
one_4pi_eps = ONE_4PI_EPS0/pme->epsilon_r;
factor = M_PI*M_PI/(pme->ewaldcoeff*pme->ewaldcoeff);
boxfactor = M_PI*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2];
esum = 0;
virxx = 0;
......@@ -444,21 +445,21 @@ pme_reciprocal_convolution(pme_t pme,
viryz = 0;
virzz = 0;
maxkx = (RealOpenMM) ((nx+1)/2);
maxky = (RealOpenMM) ((ny+1)/2);
maxkz = (RealOpenMM) ((nz+1)/2);
maxkx = (nx+1)/2;
maxky = (ny+1)/2;
maxkz = (nz+1)/2;
for (kx=0;kx<nx;kx++)
{
/* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
mx = (RealOpenMM) ((kx<maxkx) ? kx : (kx-nx));
mx = (kx<maxkx) ? kx : (kx-nx);
mhx = mx*recipBoxVectors[0][0];
bx = boxfactor*pme->bsplines_moduli[0][kx];
for (ky=0;ky<ny;ky++)
{
/* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
my = (RealOpenMM) ((ky<maxky) ? ky : (ky-ny));
my = (ky<maxky) ? ky : (ky-ny);
mhy = mx*recipBoxVectors[1][0]+my*recipBoxVectors[1][1];
by = pme->bsplines_moduli[1][ky];
......@@ -478,7 +479,7 @@ pme_reciprocal_convolution(pme_t pme,
}
/* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
mz = (RealOpenMM) ((kz<maxkz) ? kz : (kz-nz));
mz = (kz<maxkz) ? kz : (kz-nz);
mhz = mx*recipBoxVectors[2][0]+my*recipBoxVectors[2][1]+mz*recipBoxVectors[2][2];
/* Pointer to the grid cell in question */
......@@ -509,15 +510,15 @@ pme_reciprocal_convolution(pme_t pme,
}
/* The factor 0.5 is nothing special, but it is better to have it here than inside the loop :-) */
*energy = (RealOpenMM) (0.5*esum);
*energy = 0.5*esum;
}
static void
pme_grid_interpolate_force(pme_t pme,
const RealVec recipBoxVectors[3],
const vector<RealOpenMM>& charges,
vector<RealVec>& forces)
const Vec3 recipBoxVectors[3],
const vector<double>& charges,
vector<Vec3>& forces)
{
int i;
int ix,iy,iz;
......@@ -525,17 +526,17 @@ pme_grid_interpolate_force(pme_t pme,
int xindex,yindex,zindex;
int index;
int order;
RealOpenMM q;
RealOpenMM * thetax;
RealOpenMM * thetay;
RealOpenMM * thetaz;
RealOpenMM * dthetax;
RealOpenMM * dthetay;
RealOpenMM * dthetaz;
RealOpenMM tx,ty,tz;
RealOpenMM dtx,dty,dtz;
RealOpenMM fx,fy,fz;
RealOpenMM gridvalue;
double q;
double * thetax;
double * thetay;
double * thetaz;
double * dthetax;
double * dthetay;
double * dthetaz;
double tx,ty,tz;
double dtx,dty,dtz;
double fx,fy,fz;
double gridvalue;
int nx,ny,nz;
nx = pme->ngrid[0];
......@@ -617,11 +618,11 @@ pme_grid_interpolate_force(pme_t pme,
int
pme_init(pme_t * ppme,
RealOpenMM ewaldcoeff,
double ewaldcoeff,
int natoms,
const int ngrid[3],
const int ngrid[3],
int pme_order,
RealOpenMM epsilon_r)
double epsilon_r)
{
pme_t pme;
int d;
......@@ -636,8 +637,8 @@ pme_init(pme_t * ppme,
for (d=0;d<3;d++)
{
pme->ngrid[d] = ngrid[d];
pme->bsplines_theta[d] = (RealOpenMM *)malloc(sizeof(RealOpenMM)*pme_order*natoms);
pme->bsplines_dtheta[d] = (RealOpenMM *)malloc(sizeof(RealOpenMM)*pme_order*natoms);
pme->bsplines_theta[d] = (double *)malloc(sizeof(double)*pme_order*natoms);
pme->bsplines_dtheta[d] = (double *)malloc(sizeof(double)*pme_order*natoms);
}
pme->particlefraction = (rvec *)malloc(sizeof(rvec)*natoms);
......@@ -661,15 +662,15 @@ pme_init(pme_t * ppme,
int pme_exec(pme_t pme,
const vector<RealVec>& atomCoordinates,
vector<RealVec>& forces,
const vector<RealOpenMM>& charges,
const RealVec periodicBoxVectors[3],
RealOpenMM* energy)
const vector<Vec3>& atomCoordinates,
vector<Vec3>& forces,
const vector<double>& charges,
const Vec3 periodicBoxVectors[3],
double* energy)
{
/* Routine is called with coordinates in x, a box, and charges in q */
RealVec recipBoxVectors[3];
Vec3 recipBoxVectors[3];
invert_box_vectors(periodicBoxVectors, recipBoxVectors);
/* Before we can do the actual interpolation, we need to recalculate and update
......
platforms/reference/src/SimTKReference/ReferencePairIxn.cpp
View file @ a783b996
......@@ -38,13 +38,6 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */
ReferencePairIxn::ReferencePairIxn() {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferencePairIxn::ReferencePairIxn";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
......@@ -54,11 +47,4 @@ ReferencePairIxn::ReferencePairIxn() {
--------------------------------------------------------------------------------------- */
ReferencePairIxn::~ReferencePairIxn() {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferencePairIxn::~ReferencePairIxn";
// ---------------------------------------------------------------------------------------
}
platforms/reference/src/SimTKReference/ReferenceProperDihedralBond.cpp
View file @ a783b996
......@@ -50,7 +50,7 @@ ReferenceProperDihedralBond::ReferenceProperDihedralBond() : usePeriodic(false)
ReferenceProperDihedralBond::~ReferenceProperDihedralBond() {
}
void ReferenceProperDihedralBond::setPeriodic(OpenMM::RealVec* vectors) {
void ReferenceProperDihedralBond::setPeriodic(OpenMM::Vec3* vectors) {
usePeriodic = true;
boxVectors[0] = vectors[0];
boxVectors[1] = vectors[1];
......@@ -72,28 +72,13 @@ void ReferenceProperDihedralBond::setPeriodic(OpenMM::RealVec* vectors) {
--------------------------------------------------------------------------------------- */
void ReferenceProperDihedralBond::calculateBondIxn(int* atomIndices,
vector<RealVec>& atomCoordinates,
RealOpenMM* parameters,
vector<RealVec>& forces,
RealOpenMM* totalEnergy, double* energyParamDerivs) {
vector<Vec3>& atomCoordinates,
double* parameters,
vector<Vec3>& forces,
double* totalEnergy, double* energyParamDerivs) {
double deltaR[3][ReferenceForce::LastDeltaRIndex];
static const std::string methodName = "\nReferenceProperDihedralBond::calculateBondIxn";
// constants -- reduce Visual Studio warnings regarding conversions between float & double
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM oneM = -1.0;
static const int threeI = 3;
static const int LastAtomIndex = 4;
RealOpenMM deltaR[3][ReferenceForce::LastDeltaRIndex];
RealOpenMM crossProductMemory[6];
double crossProductMemory[6];
// ---------------------------------------------------------------------------------------
......@@ -114,55 +99,55 @@ void ReferenceProperDihedralBond::calculateBondIxn(int* atomIndices,
ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);
}
RealOpenMM dotDihedral;
RealOpenMM signOfAngle;
int hasREntry = 1;
double dotDihedral;
double signOfAngle;
int hasREntry = 1;
// Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
RealOpenMM* crossProduct[2];
double* crossProduct[2];
crossProduct[0] = crossProductMemory;
crossProduct[1] = crossProductMemory + 3;
// get dihedral angle
RealOpenMM dihedralAngle = getDihedralAngleBetweenThreeVectors(deltaR[0], deltaR[1], deltaR[2],
double dihedralAngle = getDihedralAngleBetweenThreeVectors(deltaR[0], deltaR[1], deltaR[2],
crossProduct, &dotDihedral, deltaR[0],
&signOfAngle, hasREntry);
// evaluate delta angle, dE/d(angle)
RealOpenMM deltaAngle = parameters[2]*dihedralAngle - parameters[1];
RealOpenMM sinDeltaAngle = SIN(deltaAngle);
RealOpenMM dEdAngle = -parameters[0]*parameters[2]*sinDeltaAngle;
RealOpenMM energy = parameters[0]*(one + COS(deltaAngle));
double deltaAngle = parameters[2]*dihedralAngle - parameters[1];
double sinDeltaAngle = SIN(deltaAngle);
double dEdAngle = -parameters[0]*parameters[2]*sinDeltaAngle;
double energy = parameters[0]*(1.0 + cos(deltaAngle));
// compute force
RealOpenMM internalF[4][3];
RealOpenMM forceFactors[4];
RealOpenMM normCross1 = DOT3(crossProduct[0], crossProduct[0]);
RealOpenMM normBC = deltaR[1][ReferenceForce::RIndex];
forceFactors[0] = (-dEdAngle*normBC)/normCross1;
double internalF[4][3];
double forceFactors[4];
double normCross1 = DOT3(crossProduct[0], crossProduct[0]);
double normBC = deltaR[1][ReferenceForce::RIndex];
forceFactors[0] = (-dEdAngle*normBC)/normCross1;
RealOpenMM normCross2 = DOT3(crossProduct[1], crossProduct[1]);
forceFactors[3] = (dEdAngle*normBC)/normCross2;
forceFactors[1] = DOT3(deltaR[0], deltaR[1]);
forceFactors[1] /= deltaR[1][ReferenceForce::R2Index];
double normCross2 = DOT3(crossProduct[1], crossProduct[1]);
forceFactors[3] = (dEdAngle*normBC)/normCross2;
forceFactors[1] = DOT3(deltaR[0], deltaR[1]);
forceFactors[1] /= deltaR[1][ReferenceForce::R2Index];
forceFactors[2] = DOT3(deltaR[2], deltaR[1]);
forceFactors[2] /= deltaR[1][ReferenceForce::R2Index];
forceFactors[2] = DOT3(deltaR[2], deltaR[1]);
forceFactors[2] /= deltaR[1][ReferenceForce::R2Index];
for (int ii = 0; ii < 3; ii++) {
internalF[0][ii] = forceFactors[0]*crossProduct[0][ii];
internalF[3][ii] = forceFactors[3]*crossProduct[1][ii];
internalF[0][ii] = forceFactors[0]*crossProduct[0][ii];
internalF[3][ii] = forceFactors[3]*crossProduct[1][ii];
RealOpenMM s = forceFactors[1]*internalF[0][ii] - forceFactors[2]*internalF[3][ii];
double s = forceFactors[1]*internalF[0][ii] - forceFactors[2]*internalF[3][ii];
internalF[1][ii] = internalF[0][ii] - s;
internalF[2][ii] = internalF[3][ii] + s;
internalF[1][ii] = internalF[0][ii] - s;
internalF[2][ii] = internalF[3][ii] + s;
}
// accumulate forces
......
platforms/reference/src/SimTKReference/ReferenceRbDihedralBond.cpp
View file @ a783b996
......@@ -50,7 +50,7 @@ ReferenceRbDihedralBond::ReferenceRbDihedralBond() : usePeriodic(false) {
ReferenceRbDihedralBond::~ReferenceRbDihedralBond() {
}
void ReferenceRbDihedralBond::setPeriodic(OpenMM::RealVec* vectors) {
void ReferenceRbDihedralBond::setPeriodic(OpenMM::Vec3* vectors) {
usePeriodic = true;
boxVectors[0] = vectors[0];
boxVectors[1] = vectors[1];
......@@ -70,32 +70,17 @@ void ReferenceRbDihedralBond::setPeriodic(OpenMM::RealVec* vectors) {
--------------------------------------------------------------------------------------- */
void ReferenceRbDihedralBond::calculateBondIxn(int* atomIndices,
vector<RealVec>& atomCoordinates,
RealOpenMM* parameters,
vector<RealVec>& forces,
RealOpenMM* totalEnergy, double* energyParamDerivs) {
static const std::string methodName = "\nReferenceRbDihedralBond::calculateBondIxn";
// constants -- reduce Visual Studio warnings regarding conversions between float & double
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM oneM = -1.0;
static const int threeI = 3;
vector<Vec3>& atomCoordinates,
double* parameters,
vector<Vec3>& forces,
double* totalEnergy, double* energyParamDerivs) {
// number of parameters
static const int numberOfParameters = 6;
static const int LastAtomIndex = 4;
RealOpenMM deltaR[3][ReferenceForce::LastDeltaRIndex];
double deltaR[3][ReferenceForce::LastDeltaRIndex];
RealOpenMM crossProductMemory[6];
double crossProductMemory[6];
// ---------------------------------------------------------------------------------------
......@@ -116,65 +101,65 @@ void ReferenceRbDihedralBond::calculateBondIxn(int* atomIndices,
ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);
}
RealOpenMM cosPhi;
RealOpenMM signOfAngle;
int hasREntry = 1;
double cosPhi;
double signOfAngle;
int hasREntry = 1;
// Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
RealOpenMM* crossProduct[2];
double* crossProduct[2];
crossProduct[0] = crossProductMemory;
crossProduct[1] = crossProductMemory + 3;
RealOpenMM dihederalAngle = getDihedralAngleBetweenThreeVectors(deltaR[0], deltaR[1], deltaR[2],
crossProduct, &cosPhi, deltaR[0],
&signOfAngle, hasREntry);
double dihederalAngle = getDihedralAngleBetweenThreeVectors(deltaR[0], deltaR[1], deltaR[2],
crossProduct, &cosPhi, deltaR[0],
&signOfAngle, hasREntry);
// Gromacs: use polymer convention
if (dihederalAngle < zero) {
if (dihederalAngle < 0.0) {
dihederalAngle += PI_M;
} else {
dihederalAngle -= PI_M;
}
cosPhi *= -one;
cosPhi *= -1.0;
// Ryckaert-Bellemans:
// V = sum over i: { C_i*cos(psi)**i }, where psi = phi - PI,
// C_i is ith RB coefficient
RealOpenMM dEdAngle = zero;
RealOpenMM energy = parameters[0];
RealOpenMM cosFactor = one;
double dEdAngle = 0.0;
double energy = parameters[0];
double cosFactor = 1.0;
for (int ii = 1; ii < numberOfParameters; ii++) {
dEdAngle -= ((RealOpenMM) ii)*parameters[ii]*cosFactor;
dEdAngle -= ii*parameters[ii]*cosFactor;
cosFactor *= cosPhi;
energy += cosFactor*parameters[ii];
}
dEdAngle *= SIN(dihederalAngle);
RealOpenMM internalF[4][3];
RealOpenMM forceFactors[4];
RealOpenMM normCross1 = DOT3(crossProduct[0], crossProduct[0]);
RealOpenMM normBC = deltaR[1][ReferenceForce::RIndex];
forceFactors[0] = (-dEdAngle*normBC)/normCross1;
double internalF[4][3];
double forceFactors[4];
double normCross1 = DOT3(crossProduct[0], crossProduct[0]);
double normBC = deltaR[1][ReferenceForce::RIndex];
forceFactors[0] = (-dEdAngle*normBC)/normCross1;
double normCross2 = DOT3(crossProduct[1], crossProduct[1]);
forceFactors[3] = (dEdAngle*normBC)/normCross2;
RealOpenMM normCross2 = DOT3(crossProduct[1], crossProduct[1]);
forceFactors[3] = (dEdAngle*normBC)/normCross2;
forceFactors[1] = DOT3(deltaR[0], deltaR[1]);
forceFactors[1] /= deltaR[1][ReferenceForce::R2Index];
forceFactors[1] = DOT3(deltaR[0], deltaR[1]);
forceFactors[1] /= deltaR[1][ReferenceForce::R2Index];
forceFactors[2] = DOT3(deltaR[2], deltaR[1]);
forceFactors[2] /= deltaR[1][ReferenceForce::R2Index];
forceFactors[2] = DOT3(deltaR[2], deltaR[1]);
forceFactors[2] /= deltaR[1][ReferenceForce::R2Index];
for (int ii = 0; ii < 3; ii++) {
internalF[0][ii] = forceFactors[0]*crossProduct[0][ii];
internalF[3][ii] = forceFactors[3]*crossProduct[1][ii];
RealOpenMM s = forceFactors[1]*internalF[0][ii] - forceFactors[2]*internalF[3][ii];
double s = forceFactors[1]*internalF[0][ii] - forceFactors[2]*internalF[3][ii];
internalF[1][ii] = internalF[0][ii] - s;
internalF[2][ii] = internalF[3][ii] + s;
......
platforms/reference/src/SimTKReference/ReferenceSETTLEAlgorithm.cpp
View file @ a783b996
......@@ -35,7 +35,7 @@ using namespace OpenMM;
using namespace std;
ReferenceSETTLEAlgorithm::ReferenceSETTLEAlgorithm(const vector<int>& atom1, const vector<int>& atom2, const vector<int>& atom3,
const vector<RealOpenMM>& distance1, const vector<RealOpenMM>& distance2, vector<RealOpenMM>& masses) :
const vector<double>& distance1, const vector<double>& distance2, vector<double>& masses) :
atom1(atom1), atom2(atom2), atom3(atom3), distance1(distance1), distance2(distance2), masses(masses) {
}
......@@ -43,7 +43,7 @@ int ReferenceSETTLEAlgorithm::getNumClusters() const {
return atom1.size();
}
void ReferenceSETTLEAlgorithm::getClusterParameters(int index, int& atom1, int& atom2, int& atom3, RealOpenMM& distance1, RealOpenMM& distance2) const {
void ReferenceSETTLEAlgorithm::getClusterParameters(int index, int& atom1, int& atom2, int& atom3, double& distance1, double& distance2) const {
atom1 = this->atom1[index];
atom2 = this->atom2[index];
atom3 = this->atom3[index];
......@@ -51,130 +51,130 @@ void ReferenceSETTLEAlgorithm::getClusterParameters(int index, int& atom1, int&
distance2 = this->distance2[index];
}
void ReferenceSETTLEAlgorithm::apply(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) {
void ReferenceSETTLEAlgorithm::apply(vector<OpenMM::Vec3>& atomCoordinates, vector<OpenMM::Vec3>& atomCoordinatesP, vector<double>& inverseMasses, double tolerance) {
for (int index = 0; index < (int) atom1.size(); ++index) {
RealVec apos0 = atomCoordinates[atom1[index]];
RealVec xp0 = atomCoordinatesP[atom1[index]]-apos0;
RealVec apos1 = atomCoordinates[atom2[index]];
RealVec xp1 = atomCoordinatesP[atom2[index]]-apos1;
RealVec apos2 = atomCoordinates[atom3[index]];
RealVec xp2 = atomCoordinatesP[atom3[index]]-apos2;
RealOpenMM m0 = masses[atom1[index]];
RealOpenMM m1 = masses[atom2[index]];
RealOpenMM m2 = masses[atom3[index]];
Vec3 apos0 = atomCoordinates[atom1[index]];
Vec3 xp0 = atomCoordinatesP[atom1[index]]-apos0;
Vec3 apos1 = atomCoordinates[atom2[index]];
Vec3 xp1 = atomCoordinatesP[atom2[index]]-apos1;
Vec3 apos2 = atomCoordinates[atom3[index]];
Vec3 xp2 = atomCoordinatesP[atom3[index]]-apos2;
double m0 = masses[atom1[index]];
double m1 = masses[atom2[index]];
double m2 = masses[atom3[index]];
// Apply the SETTLE algorithm.
RealOpenMM xb0 = apos1[0]-apos0[0];
RealOpenMM yb0 = apos1[1]-apos0[1];
RealOpenMM zb0 = apos1[2]-apos0[2];
RealOpenMM xc0 = apos2[0]-apos0[0];
RealOpenMM yc0 = apos2[1]-apos0[1];
RealOpenMM zc0 = apos2[2]-apos0[2];
RealOpenMM invTotalMass = 1/(m0+m1+m2);
RealOpenMM xcom = (xp0[0]*m0 + (xb0+xp1[0])*m1 + (xc0+xp2[0])*m2) * invTotalMass;
RealOpenMM ycom = (xp0[1]*m0 + (yb0+xp1[1])*m1 + (yc0+xp2[1])*m2) * invTotalMass;
RealOpenMM zcom = (xp0[2]*m0 + (zb0+xp1[2])*m1 + (zc0+xp2[2])*m2) * invTotalMass;
RealOpenMM xa1 = xp0[0] - xcom;
RealOpenMM ya1 = xp0[1] - ycom;
RealOpenMM za1 = xp0[2] - zcom;
RealOpenMM xb1 = xb0 + xp1[0] - xcom;
RealOpenMM yb1 = yb0 + xp1[1] - ycom;
RealOpenMM zb1 = zb0 + xp1[2] - zcom;
RealOpenMM xc1 = xc0 + xp2[0] - xcom;
RealOpenMM yc1 = yc0 + xp2[1] - ycom;
RealOpenMM zc1 = zc0 + xp2[2] - zcom;
RealOpenMM xaksZd = yb0*zc0 - zb0*yc0;
RealOpenMM yaksZd = zb0*xc0 - xb0*zc0;
RealOpenMM zaksZd = xb0*yc0 - yb0*xc0;
RealOpenMM xaksXd = ya1*zaksZd - za1*yaksZd;
RealOpenMM yaksXd = za1*xaksZd - xa1*zaksZd;
RealOpenMM zaksXd = xa1*yaksZd - ya1*xaksZd;
RealOpenMM xaksYd = yaksZd*zaksXd - zaksZd*yaksXd;
RealOpenMM yaksYd = zaksZd*xaksXd - xaksZd*zaksXd;
RealOpenMM zaksYd = xaksZd*yaksXd - yaksZd*xaksXd;
RealOpenMM axlng = sqrt(xaksXd*xaksXd + yaksXd*yaksXd + zaksXd*zaksXd);
RealOpenMM aylng = sqrt(xaksYd*xaksYd + yaksYd*yaksYd + zaksYd*zaksYd);
RealOpenMM azlng = sqrt(xaksZd*xaksZd + yaksZd*yaksZd + zaksZd*zaksZd);
RealOpenMM trns11 = xaksXd / axlng;
RealOpenMM trns21 = yaksXd / axlng;
RealOpenMM trns31 = zaksXd / axlng;
RealOpenMM trns12 = xaksYd / aylng;
RealOpenMM trns22 = yaksYd / aylng;
RealOpenMM trns32 = zaksYd / aylng;
RealOpenMM trns13 = xaksZd / azlng;
RealOpenMM trns23 = yaksZd / azlng;
RealOpenMM trns33 = zaksZd / azlng;
RealOpenMM xb0d = trns11*xb0 + trns21*yb0 + trns31*zb0;
RealOpenMM yb0d = trns12*xb0 + trns22*yb0 + trns32*zb0;
RealOpenMM xc0d = trns11*xc0 + trns21*yc0 + trns31*zc0;
RealOpenMM yc0d = trns12*xc0 + trns22*yc0 + trns32*zc0;
RealOpenMM za1d = trns13*xa1 + trns23*ya1 + trns33*za1;
RealOpenMM xb1d = trns11*xb1 + trns21*yb1 + trns31*zb1;
RealOpenMM yb1d = trns12*xb1 + trns22*yb1 + trns32*zb1;
RealOpenMM zb1d = trns13*xb1 + trns23*yb1 + trns33*zb1;
RealOpenMM xc1d = trns11*xc1 + trns21*yc1 + trns31*zc1;
RealOpenMM yc1d = trns12*xc1 + trns22*yc1 + trns32*zc1;
RealOpenMM zc1d = trns13*xc1 + trns23*yc1 + trns33*zc1;
double xb0 = apos1[0]-apos0[0];
double yb0 = apos1[1]-apos0[1];
double zb0 = apos1[2]-apos0[2];
double xc0 = apos2[0]-apos0[0];
double yc0 = apos2[1]-apos0[1];
double zc0 = apos2[2]-apos0[2];
double invTotalMass = 1/(m0+m1+m2);
double xcom = (xp0[0]*m0 + (xb0+xp1[0])*m1 + (xc0+xp2[0])*m2) * invTotalMass;
double ycom = (xp0[1]*m0 + (yb0+xp1[1])*m1 + (yc0+xp2[1])*m2) * invTotalMass;
double zcom = (xp0[2]*m0 + (zb0+xp1[2])*m1 + (zc0+xp2[2])*m2) * invTotalMass;
double xa1 = xp0[0] - xcom;
double ya1 = xp0[1] - ycom;
double za1 = xp0[2] - zcom;
double xb1 = xb0 + xp1[0] - xcom;
double yb1 = yb0 + xp1[1] - ycom;
double zb1 = zb0 + xp1[2] - zcom;
double xc1 = xc0 + xp2[0] - xcom;
double yc1 = yc0 + xp2[1] - ycom;
double zc1 = zc0 + xp2[2] - zcom;
double xaksZd = yb0*zc0 - zb0*yc0;
double yaksZd = zb0*xc0 - xb0*zc0;
double zaksZd = xb0*yc0 - yb0*xc0;
double xaksXd = ya1*zaksZd - za1*yaksZd;
double yaksXd = za1*xaksZd - xa1*zaksZd;
double zaksXd = xa1*yaksZd - ya1*xaksZd;
double xaksYd = yaksZd*zaksXd - zaksZd*yaksXd;
double yaksYd = zaksZd*xaksXd - xaksZd*zaksXd;
double zaksYd = xaksZd*yaksXd - yaksZd*xaksXd;
double axlng = sqrt(xaksXd*xaksXd + yaksXd*yaksXd + zaksXd*zaksXd);
double aylng = sqrt(xaksYd*xaksYd + yaksYd*yaksYd + zaksYd*zaksYd);
double azlng = sqrt(xaksZd*xaksZd + yaksZd*yaksZd + zaksZd*zaksZd);
double trns11 = xaksXd / axlng;
double trns21 = yaksXd / axlng;
double trns31 = zaksXd / axlng;
double trns12 = xaksYd / aylng;
double trns22 = yaksYd / aylng;
double trns32 = zaksYd / aylng;
double trns13 = xaksZd / azlng;
double trns23 = yaksZd / azlng;
double trns33 = zaksZd / azlng;
double xb0d = trns11*xb0 + trns21*yb0 + trns31*zb0;
double yb0d = trns12*xb0 + trns22*yb0 + trns32*zb0;
double xc0d = trns11*xc0 + trns21*yc0 + trns31*zc0;
double yc0d = trns12*xc0 + trns22*yc0 + trns32*zc0;
double za1d = trns13*xa1 + trns23*ya1 + trns33*za1;
double xb1d = trns11*xb1 + trns21*yb1 + trns31*zb1;
double yb1d = trns12*xb1 + trns22*yb1 + trns32*zb1;
double zb1d = trns13*xb1 + trns23*yb1 + trns33*zb1;
double xc1d = trns11*xc1 + trns21*yc1 + trns31*zc1;
double yc1d = trns12*xc1 + trns22*yc1 + trns32*zc1;
double zc1d = trns13*xc1 + trns23*yc1 + trns33*zc1;
// --- Step2 A2' ---
RealOpenMM rc = 0.5*distance2[index];
RealOpenMM rb = sqrt(distance1[index]*distance1[index]-rc*rc);
RealOpenMM ra = rb*(m1+m2)*invTotalMass;
double rc = 0.5*distance2[index];
double rb = sqrt(distance1[index]*distance1[index]-rc*rc);
double ra = rb*(m1+m2)*invTotalMass;
rb -= ra;
RealOpenMM sinphi = za1d / ra;
RealOpenMM cosphi = sqrt(1 - sinphi*sinphi);
RealOpenMM sinpsi = (zb1d - zc1d) / (2*rc*cosphi);
RealOpenMM cospsi = sqrt(1 - sinpsi*sinpsi);
RealOpenMM ya2d = ra*cosphi;
RealOpenMM xb2d = - rc*cospsi;
RealOpenMM yb2d = - rb*cosphi - rc*sinpsi*sinphi;
RealOpenMM yc2d = - rb*cosphi + rc*sinpsi*sinphi;
RealOpenMM xb2d2 = xb2d*xb2d;
RealOpenMM hh2 = 4.0f*xb2d2 + (yb2d-yc2d)*(yb2d-yc2d) + (zb1d-zc1d)*(zb1d-zc1d);
RealOpenMM deltx = 2.0f*xb2d + sqrt(4.0f*xb2d2 - hh2 + distance2[index]*distance2[index]);
double sinphi = za1d / ra;
double cosphi = sqrt(1 - sinphi*sinphi);
double sinpsi = (zb1d - zc1d) / (2*rc*cosphi);
double cospsi = sqrt(1 - sinpsi*sinpsi);
double ya2d = ra*cosphi;
double xb2d = - rc*cospsi;
double yb2d = - rb*cosphi - rc*sinpsi*sinphi;
double yc2d = - rb*cosphi + rc*sinpsi*sinphi;
double xb2d2 = xb2d*xb2d;
double hh2 = 4.0f*xb2d2 + (yb2d-yc2d)*(yb2d-yc2d) + (zb1d-zc1d)*(zb1d-zc1d);
double deltx = 2.0f*xb2d + sqrt(4.0f*xb2d2 - hh2 + distance2[index]*distance2[index]);
xb2d -= deltx*0.5;
// --- Step3 al,be,ga ---
RealOpenMM alpha = (xb2d*(xb0d-xc0d) + yb0d*yb2d + yc0d*yc2d);
RealOpenMM beta = (xb2d*(yc0d-yb0d) + xb0d*yb2d + xc0d*yc2d);
RealOpenMM gamma = xb0d*yb1d - xb1d*yb0d + xc0d*yc1d - xc1d*yc0d;
double alpha = (xb2d*(xb0d-xc0d) + yb0d*yb2d + yc0d*yc2d);
double beta = (xb2d*(yc0d-yb0d) + xb0d*yb2d + xc0d*yc2d);
double gamma = xb0d*yb1d - xb1d*yb0d + xc0d*yc1d - xc1d*yc0d;
RealOpenMM al2be2 = alpha*alpha + beta*beta;
RealOpenMM sintheta = (alpha*gamma - beta*sqrt(al2be2 - gamma*gamma)) / al2be2;
double al2be2 = alpha*alpha + beta*beta;
double sintheta = (alpha*gamma - beta*sqrt(al2be2 - gamma*gamma)) / al2be2;
// --- Step4 A3' ---
RealOpenMM costheta = sqrt(1 - sintheta*sintheta);
RealOpenMM xa3d = - ya2d*sintheta;
RealOpenMM ya3d = ya2d*costheta;
RealOpenMM za3d = za1d;
RealOpenMM xb3d = xb2d*costheta - yb2d*sintheta;
RealOpenMM yb3d = xb2d*sintheta + yb2d*costheta;
RealOpenMM zb3d = zb1d;
RealOpenMM xc3d = - xb2d*costheta - yc2d*sintheta;
RealOpenMM yc3d = - xb2d*sintheta + yc2d*costheta;
RealOpenMM zc3d = zc1d;
double costheta = sqrt(1 - sintheta*sintheta);
double xa3d = - ya2d*sintheta;
double ya3d = ya2d*costheta;
double za3d = za1d;
double xb3d = xb2d*costheta - yb2d*sintheta;
double yb3d = xb2d*sintheta + yb2d*costheta;
double zb3d = zb1d;
double xc3d = - xb2d*costheta - yc2d*sintheta;
double yc3d = - xb2d*sintheta + yc2d*costheta;
double zc3d = zc1d;
// --- Step5 A3 ---
RealOpenMM xa3 = trns11*xa3d + trns12*ya3d + trns13*za3d;
RealOpenMM ya3 = trns21*xa3d + trns22*ya3d + trns23*za3d;
RealOpenMM za3 = trns31*xa3d + trns32*ya3d + trns33*za3d;
RealOpenMM xb3 = trns11*xb3d + trns12*yb3d + trns13*zb3d;
RealOpenMM yb3 = trns21*xb3d + trns22*yb3d + trns23*zb3d;
RealOpenMM zb3 = trns31*xb3d + trns32*yb3d + trns33*zb3d;
RealOpenMM xc3 = trns11*xc3d + trns12*yc3d + trns13*zc3d;
RealOpenMM yc3 = trns21*xc3d + trns22*yc3d + trns23*zc3d;
RealOpenMM zc3 = trns31*xc3d + trns32*yc3d + trns33*zc3d;
double xa3 = trns11*xa3d + trns12*ya3d + trns13*za3d;
double ya3 = trns21*xa3d + trns22*ya3d + trns23*za3d;
double za3 = trns31*xa3d + trns32*ya3d + trns33*za3d;
double xb3 = trns11*xb3d + trns12*yb3d + trns13*zb3d;
double yb3 = trns21*xb3d + trns22*yb3d + trns23*zb3d;
double zb3 = trns31*xb3d + trns32*yb3d + trns33*zb3d;
double xc3 = trns11*xc3d + trns12*yc3d + trns13*zc3d;
double yc3 = trns21*xc3d + trns22*yc3d + trns23*zc3d;
double zc3 = trns31*xc3d + trns32*yc3d + trns33*zc3d;
xp0[0] = xcom + xa3;
xp0[1] = ycom + ya3;
......@@ -194,46 +194,46 @@ void ReferenceSETTLEAlgorithm::apply(vector<OpenMM::RealVec>& atomCoordinates, v
}
}
void ReferenceSETTLEAlgorithm::applyToVelocities(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& velocities, vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) {
void ReferenceSETTLEAlgorithm::applyToVelocities(vector<OpenMM::Vec3>& atomCoordinates, vector<OpenMM::Vec3>& velocities, vector<double>& inverseMasses, double tolerance) {
for (int index = 0; index < (int) atom1.size(); ++index) {
RealVec apos0 = atomCoordinates[atom1[index]];
RealVec apos1 = atomCoordinates[atom2[index]];
RealVec apos2 = atomCoordinates[atom3[index]];
RealVec v0 = velocities[atom1[index]];
RealVec v1 = velocities[atom2[index]];
RealVec v2 = velocities[atom3[index]];
Vec3 apos0 = atomCoordinates[atom1[index]];
Vec3 apos1 = atomCoordinates[atom2[index]];
Vec3 apos2 = atomCoordinates[atom3[index]];
Vec3 v0 = velocities[atom1[index]];
Vec3 v1 = velocities[atom2[index]];
Vec3 v2 = velocities[atom3[index]];
// Compute intermediate quantities: the atom masses, the bond directions, the relative velocities,
// and the angle cosines and sines.
RealOpenMM mA = masses[atom1[index]];
RealOpenMM mB = masses[atom2[index]];
RealOpenMM mC = masses[atom3[index]];
RealVec eAB = apos1-apos0;
RealVec eBC = apos2-apos1;
RealVec eCA = apos0-apos2;
double mA = masses[atom1[index]];
double mB = masses[atom2[index]];
double mC = masses[atom3[index]];
Vec3 eAB = apos1-apos0;
Vec3 eBC = apos2-apos1;
Vec3 eCA = apos0-apos2;
eAB /= sqrt(eAB[0]*eAB[0] + eAB[1]*eAB[1] + eAB[2]*eAB[2]);
eBC /= sqrt(eBC[0]*eBC[0] + eBC[1]*eBC[1] + eBC[2]*eBC[2]);
eCA /= sqrt(eCA[0]*eCA[0] + eCA[1]*eCA[1] + eCA[2]*eCA[2]);
RealOpenMM vAB = (v1[0]-v0[0])*eAB[0] + (v1[1]-v0[1])*eAB[1] + (v1[2]-v0[2])*eAB[2];
RealOpenMM vBC = (v2[0]-v1[0])*eBC[0] + (v2[1]-v1[1])*eBC[1] + (v2[2]-v1[2])*eBC[2];
RealOpenMM vCA = (v0[0]-v2[0])*eCA[0] + (v0[1]-v2[1])*eCA[1] + (v0[2]-v2[2])*eCA[2];
RealOpenMM cA = -(eAB[0]*eCA[0] + eAB[1]*eCA[1] + eAB[2]*eCA[2]);
RealOpenMM cB = -(eAB[0]*eBC[0] + eAB[1]*eBC[1] + eAB[2]*eBC[2]);
RealOpenMM cC = -(eBC[0]*eCA[0] + eBC[1]*eCA[1] + eBC[2]*eCA[2]);
RealOpenMM s2A = 1-cA*cA;
RealOpenMM s2B = 1-cB*cB;
RealOpenMM s2C = 1-cC*cC;
double vAB = (v1[0]-v0[0])*eAB[0] + (v1[1]-v0[1])*eAB[1] + (v1[2]-v0[2])*eAB[2];
double vBC = (v2[0]-v1[0])*eBC[0] + (v2[1]-v1[1])*eBC[1] + (v2[2]-v1[2])*eBC[2];
double vCA = (v0[0]-v2[0])*eCA[0] + (v0[1]-v2[1])*eCA[1] + (v0[2]-v2[2])*eCA[2];
double cA = -(eAB[0]*eCA[0] + eAB[1]*eCA[1] + eAB[2]*eCA[2]);
double cB = -(eAB[0]*eBC[0] + eAB[1]*eBC[1] + eAB[2]*eBC[2]);
double cC = -(eBC[0]*eCA[0] + eBC[1]*eCA[1] + eBC[2]*eCA[2]);
double s2A = 1-cA*cA;
double s2B = 1-cB*cB;
double s2C = 1-cC*cC;
// Solve the equations. These are different from those in the SETTLE paper (JCC 13(8), pp. 952-962, 1992), because
// in going from equations B1 to B2, they make the assumption that mB=mC (but don't bother to mention they're
// making that assumption). We allow all three atoms to have different masses.
RealOpenMM mABCinv = 1/(mA*mB*mC);
RealOpenMM denom = (((s2A*mB+s2B*mA)*mC+(s2A*mB*mB+2*(cA*cB*cC+1)*mA*mB+s2B*mA*mA))*mC+s2C*mA*mB*(mA+mB))*mABCinv;
RealOpenMM tab = ((cB*cC*mA-cA*mB-cA*mC)*vCA + (cA*cC*mB-cB*mC-cB*mA)*vBC + (s2C*mA*mA*mB*mB*mABCinv+(mA+mB+mC))*vAB)/denom;
RealOpenMM tbc = ((cA*cB*mC-cC*mB-cC*mA)*vCA + (s2A*mB*mB*mC*mC*mABCinv+(mA+mB+mC))*vBC + (cA*cC*mB-cB*mA-cB*mC)*vAB)/denom;
RealOpenMM tca = ((s2B*mA*mA*mC*mC*mABCinv+(mA+mB+mC))*vCA + (cA*cB*mC-cC*mB-cC*mA)*vBC + (cB*cC*mA-cA*mB-cA*mC)*vAB)/denom;
double mABCinv = 1/(mA*mB*mC);
double denom = (((s2A*mB+s2B*mA)*mC+(s2A*mB*mB+2*(cA*cB*cC+1)*mA*mB+s2B*mA*mA))*mC+s2C*mA*mB*(mA+mB))*mABCinv;
double tab = ((cB*cC*mA-cA*mB-cA*mC)*vCA + (cA*cC*mB-cB*mC-cB*mA)*vBC + (s2C*mA*mA*mB*mB*mABCinv+(mA+mB+mC))*vAB)/denom;
double tbc = ((cA*cB*mC-cC*mB-cC*mA)*vCA + (s2A*mB*mB*mC*mC*mABCinv+(mA+mB+mC))*vBC + (cA*cC*mB-cB*mA-cB*mC)*vAB)/denom;
double tca = ((s2B*mA*mA*mC*mC*mABCinv+(mA+mB+mC))*vCA + (cA*cB*mC-cC*mB-cC*mA)*vBC + (cB*cC*mA-cA*mB-cA*mC)*vAB)/denom;
v0 += (eAB*tab - eCA*tca)*inverseMasses[atom1[index]];
v1 += (eBC*tbc - eAB*tab)*inverseMasses[atom2[index]];
v2 += (eCA*tca - eBC*tbc)*inverseMasses[atom3[index]];
......
platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp
View file @ a783b996
......@@ -47,8 +47,8 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */
ReferenceStochasticDynamics::ReferenceStochasticDynamics(int numberOfAtoms,
RealOpenMM deltaT, RealOpenMM friction,
RealOpenMM temperature) :
double deltaT, double friction,
double temperature) :
ReferenceDynamics(numberOfAtoms, deltaT, temperature), friction(friction) {
xPrime.resize(numberOfAtoms);
inverseMasses.resize(numberOfAtoms);
......@@ -61,13 +61,6 @@ ReferenceStochasticDynamics::ReferenceStochasticDynamics(int numberOfAtoms,
--------------------------------------------------------------------------------------- */
ReferenceStochasticDynamics::~ReferenceStochasticDynamics() {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceStochasticDynamics::~ReferenceStochasticDynamics";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
......@@ -76,7 +69,7 @@ ReferenceStochasticDynamics::~ReferenceStochasticDynamics() {
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceStochasticDynamics::getFriction() const {
double ReferenceStochasticDynamics::getFriction() const {
return friction;
}
......@@ -93,29 +86,22 @@ RealOpenMM ReferenceStochasticDynamics::getFriction() const {
--------------------------------------------------------------------------------------- */
void ReferenceStochasticDynamics::updatePart1(int numberOfAtoms, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& inverseMasses,
vector<RealVec>& xPrime) {
// ---------------------------------------------------------------------------------------
//static const char* methodName = "\nReferenceStochasticDynamics::updatePart1";
// ---------------------------------------------------------------------------------------
void ReferenceStochasticDynamics::updatePart1(int numberOfAtoms, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities,
vector<Vec3>& forces, vector<double>& inverseMasses,
vector<Vec3>& xPrime) {
// perform first update
RealOpenMM dt = getDeltaT();
RealOpenMM friction = getFriction();
const RealOpenMM vscale = EXP(-dt*friction);
const RealOpenMM fscale = (friction == 0 ? dt : (1-vscale)/friction);
const RealOpenMM kT = BOLTZ*getTemperature();
const RealOpenMM noisescale = SQRT(kT*(1-vscale*vscale));
double dt = getDeltaT();
double friction = getFriction();
const double vscale = exp(-dt*friction);
const double fscale = (friction == 0 ? dt : (1-vscale)/friction);
const double kT = BOLTZ*getTemperature();
const double noisescale = sqrt(kT*(1-vscale*vscale));
for (int ii = 0; ii < numberOfAtoms; ii++) {
if (inverseMasses[ii] != 0.0) {
RealOpenMM sqrtInvMass = SQRT(inverseMasses[ii]);
double sqrtInvMass = sqrt(inverseMasses[ii]);
for (int jj = 0; jj < 3; jj++) {
velocities[ii][jj] = vscale*velocities[ii][jj] + fscale*inverseMasses[ii]*forces[ii][jj] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
......@@ -135,17 +121,10 @@ void ReferenceStochasticDynamics::updatePart1(int numberOfAtoms, vector<RealVec>
--------------------------------------------------------------------------------------- */
void ReferenceStochasticDynamics::updatePart2(int numberOfAtoms, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& inverseMasses,
vector<RealVec>& xPrime) {
// ---------------------------------------------------------------------------------------
//static const char* methodName = "\nReferenceStochasticDynamics::updatePart2";
// ---------------------------------------------------------------------------------------
void ReferenceStochasticDynamics::updatePart2(int numberOfAtoms, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities,
vector<Vec3>& forces, vector<double>& inverseMasses,
vector<Vec3>& xPrime) {
// perform second update
for (int ii = 0; ii < numberOfAtoms; ii++) {
......@@ -154,10 +133,10 @@ void ReferenceStochasticDynamics::updatePart2(int numberOfAtoms, vector<RealVec>
}
}
void ReferenceStochasticDynamics::updatePart3(int numberOfAtoms, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities, vector<RealOpenMM>& inverseMasses,
vector<RealVec>& xPrime) {
RealOpenMM invStepSize = 1.0/getDeltaT();
void ReferenceStochasticDynamics::updatePart3(int numberOfAtoms, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities, vector<double>& inverseMasses,
vector<Vec3>& xPrime) {
double invStepSize = 1.0/getDeltaT();
for (int i = 0; i < numberOfAtoms; ++i)
if (inverseMasses[i] != 0) {
velocities[i] = (xPrime[i]-atomCoordinates[i])*invStepSize;
......@@ -178,18 +157,8 @@ void ReferenceStochasticDynamics::updatePart3(int numberOfAtoms, vector<RealVec>
--------------------------------------------------------------------------------------- */
void ReferenceStochasticDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM tolerance) {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceStochasticDynamics::update";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
void ReferenceStochasticDynamics::update(const OpenMM::System& system, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities, vector<Vec3>& forces, vector<double>& masses, double tolerance) {
// first-time-through initialization
int numberOfAtoms = system.getNumParticles();
......@@ -197,10 +166,10 @@ void ReferenceStochasticDynamics::update(const OpenMM::System& system, vector<Re
// invert masses
for (int ii = 0; ii < numberOfAtoms; ii++) {
if (masses[ii] == zero)
inverseMasses[ii] = zero;
if (masses[ii] == 0.0)
inverseMasses[ii] = 0.0;
else
inverseMasses[ii] = one/masses[ii];
inverseMasses[ii] = 1.0/masses[ii];
}
}
......
platforms/reference/src/SimTKReference/ReferenceVariableStochasticDynamics.cpp
View file @ a783b996
......@@ -49,8 +49,8 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */
ReferenceVariableStochasticDynamics::ReferenceVariableStochasticDynamics(int numberOfAtoms,
RealOpenMM friction, RealOpenMM temperature,
RealOpenMM accuracy) :
double friction, double temperature,
double accuracy) :
ReferenceDynamics(numberOfAtoms, 0.0f, temperature), friction(friction), _accuracy(accuracy) {
xPrime.resize(numberOfAtoms);
inverseMasses.resize(numberOfAtoms);
......@@ -63,13 +63,6 @@ ReferenceVariableStochasticDynamics::ReferenceVariableStochasticDynamics(int num
--------------------------------------------------------------------------------------- */
ReferenceVariableStochasticDynamics::~ReferenceVariableStochasticDynamics() {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceVariableStochasticDynamics::~ReferenceVariableStochasticDynamics";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
......@@ -80,7 +73,7 @@ ReferenceVariableStochasticDynamics::~ReferenceVariableStochasticDynamics() {
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceVariableStochasticDynamics::getAccuracy() const {
double ReferenceVariableStochasticDynamics::getAccuracy() const {
return _accuracy;
}
......@@ -90,7 +83,7 @@ RealOpenMM ReferenceVariableStochasticDynamics::getAccuracy() const {
--------------------------------------------------------------------------------------- */
void ReferenceVariableStochasticDynamics::setAccuracy(RealOpenMM accuracy) {
void ReferenceVariableStochasticDynamics::setAccuracy(double accuracy) {
_accuracy = accuracy;
}
......@@ -100,7 +93,7 @@ void ReferenceVariableStochasticDynamics::setAccuracy(RealOpenMM accuracy) {
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceVariableStochasticDynamics::getFriction() const {
double ReferenceVariableStochasticDynamics::getFriction() const {
return friction;
}
......@@ -119,18 +112,10 @@ RealOpenMM ReferenceVariableStochasticDynamics::getFriction() const {
--------------------------------------------------------------------------------------- */
void ReferenceVariableStochasticDynamics::updatePart1(int numberOfAtoms, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& masses, vector<RealOpenMM>& inverseMasses,
vector<RealVec>& xPrime, RealOpenMM maxStepSize) {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceVariableStochasticDynamics::updatePart1";
// ---------------------------------------------------------------------------------------
void ReferenceVariableStochasticDynamics::updatePart1(int numberOfAtoms, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities,
vector<Vec3>& forces, vector<double>& masses, vector<double>& inverseMasses,
vector<Vec3>& xPrime, double maxStepSize) {
// first-time-through initialization
if (getTimeStep() == 0) {
......@@ -145,15 +130,15 @@ void ReferenceVariableStochasticDynamics::updatePart1(int numberOfAtoms, vector<
}
// Select the step size to use
RealOpenMM error = 0;
double error = 0;
for (int i = 0; i < numberOfAtoms; ++i) {
for (int j = 0; j < 3; ++j) {
RealOpenMM xerror = inverseMasses[i]*forces[i][j];
double xerror = inverseMasses[i]*forces[i][j];
error += xerror*xerror;
}
}
error = SQRT(error/(numberOfAtoms*3));
RealOpenMM newStepSize = SQRT(getAccuracy()/error);
error = sqrt(error/(numberOfAtoms*3));
double newStepSize = sqrt(getAccuracy()/error);
if (getDeltaT() > 0.0f)
newStepSize = std::min(newStepSize, getDeltaT()*2.0f); // For safety, limit how quickly dt can increase.
if (newStepSize > getDeltaT() && newStepSize < 1.2f*getDeltaT())
......@@ -164,16 +149,16 @@ void ReferenceVariableStochasticDynamics::updatePart1(int numberOfAtoms, vector<
// perform first update
RealOpenMM dt = getDeltaT();
RealOpenMM friction = getFriction();
const RealOpenMM vscale = EXP(-dt*friction);
const RealOpenMM fscale = (friction == 0 ? dt : (1-vscale)/friction);
const RealOpenMM kT = BOLTZ*getTemperature();
const RealOpenMM noisescale = SQRT(kT*(1-vscale*vscale));
double dt = getDeltaT();
double friction = getFriction();
const double vscale = exp(-dt*friction);
const double fscale = (friction == 0 ? dt : (1-vscale)/friction);
const double kT = BOLTZ*getTemperature();
const double noisescale = sqrt(kT*(1-vscale*vscale));
for (int ii = 0; ii < numberOfAtoms; ii++) {
if (masses[ii] != 0) {
RealOpenMM sqrtInvMass = SQRT(inverseMasses[ii]);
double sqrtInvMass = sqrt(inverseMasses[ii]);
for (int jj = 0; jj < 3; jj++) {
velocities[ii][jj] = vscale*velocities[ii][jj] + fscale*inverseMasses[ii]*forces[ii][jj] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
}
......@@ -194,17 +179,10 @@ void ReferenceVariableStochasticDynamics::updatePart1(int numberOfAtoms, vector<
--------------------------------------------------------------------------------------- */
void ReferenceVariableStochasticDynamics::updatePart2(int numberOfAtoms, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& inverseMasses,
vector<RealVec>& xPrime) {
// ---------------------------------------------------------------------------------------
//static const char* methodName = "\nReferenceVariableStochasticDynamics::updatePart2";
// ---------------------------------------------------------------------------------------
void ReferenceVariableStochasticDynamics::updatePart2(int numberOfAtoms, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities,
vector<Vec3>& forces, vector<double>& inverseMasses,
vector<Vec3>& xPrime) {
// perform second update
for (int ii = 0; ii < numberOfAtoms; ii++) {
......@@ -228,15 +206,9 @@ void ReferenceVariableStochasticDynamics::updatePart2(int numberOfAtoms, vector<
--------------------------------------------------------------------------------------- */
void ReferenceVariableStochasticDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM maxStepSize, RealOpenMM tolerance) {
// ---------------------------------------------------------------------------------------
//static const char* methodName = "\nReferenceVariableStochasticDynamics::update";
// ---------------------------------------------------------------------------------------
void ReferenceVariableStochasticDynamics::update(const OpenMM::System& system, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities,
vector<Vec3>& forces, vector<double>& masses, double maxStepSize, double tolerance) {
// 1st update
......@@ -253,7 +225,7 @@ void ReferenceVariableStochasticDynamics::update(const OpenMM::System& system, v
// copy xPrime -> atomCoordinates
RealOpenMM invStepSize = 1.0/getDeltaT();
double invStepSize = 1.0/getDeltaT();
for (int ii = 0; ii < numberOfAtoms; ii++) {
if (masses[ii] != 0.0) {
velocities[ii] = (xPrime[ii]-atomCoordinates[ii])*invStepSize;
......
platforms/reference/src/SimTKReference/ReferenceVariableVerletDynamics.cpp
View file @ a783b996
......@@ -43,7 +43,7 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */
ReferenceVariableVerletDynamics::ReferenceVariableVerletDynamics(int numberOfAtoms, RealOpenMM accuracy) :
ReferenceVariableVerletDynamics::ReferenceVariableVerletDynamics(int numberOfAtoms, double accuracy) :
ReferenceDynamics(numberOfAtoms, 0.0f, 0.0f), _accuracy(accuracy) {
xPrime.resize(numberOfAtoms);
inverseMasses.resize(numberOfAtoms);
......@@ -56,13 +56,6 @@ ReferenceVariableVerletDynamics::ReferenceVariableVerletDynamics(int numberOfAto
--------------------------------------------------------------------------------------- */
ReferenceVariableVerletDynamics::~ReferenceVariableVerletDynamics() {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceVariableVerletDynamics::~ReferenceVariableVerletDynamics";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
......@@ -73,7 +66,7 @@ ReferenceVariableVerletDynamics::~ReferenceVariableVerletDynamics() {
--------------------------------------------------------------------------------------- */
RealOpenMM ReferenceVariableVerletDynamics::getAccuracy() const {
double ReferenceVariableVerletDynamics::getAccuracy() const {
return _accuracy;
}
......@@ -83,7 +76,7 @@ RealOpenMM ReferenceVariableVerletDynamics::getAccuracy() const {
--------------------------------------------------------------------------------------- */
void ReferenceVariableVerletDynamics::setAccuracy(RealOpenMM accuracy) {
void ReferenceVariableVerletDynamics::setAccuracy(double accuracy) {
_accuracy = accuracy;
}
......@@ -101,19 +94,9 @@ void ReferenceVariableVerletDynamics::setAccuracy(RealOpenMM accuracy) {
--------------------------------------------------------------------------------------- */
void ReferenceVariableVerletDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM maxStepSize, RealOpenMM tolerance) {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceVariableVerletDynamics::update";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
void ReferenceVariableVerletDynamics::update(const OpenMM::System& system, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities,
vector<Vec3>& forces, vector<double>& masses, double maxStepSize, double tolerance) {
// first-time-through initialization
int numberOfAtoms = system.getNumParticles();
......@@ -121,34 +104,34 @@ void ReferenceVariableVerletDynamics::update(const OpenMM::System& system, vecto
// invert masses
for (int ii = 0; ii < numberOfAtoms; ii++) {
if (masses[ii] == zero)
inverseMasses[ii] = zero;
if (masses[ii] == 0.0)
inverseMasses[ii] = 0.0;
else
inverseMasses[ii] = one/masses[ii];
inverseMasses[ii] = 1.0/masses[ii];
}
}
RealOpenMM error = zero;
double error = 0.0;
for (int i = 0; i < numberOfAtoms; ++i) {
for (int j = 0; j < 3; ++j) {
RealOpenMM xerror = inverseMasses[i]*forces[i][j];
double xerror = inverseMasses[i]*forces[i][j];
error += xerror*xerror;
}
}
error = SQRT(error/(numberOfAtoms*3));
RealOpenMM newStepSize = SQRT(getAccuracy()/error);
error = sqrt(error/(numberOfAtoms*3));
double newStepSize = sqrt(getAccuracy()/error);
if (getDeltaT() > 0.0f)
newStepSize = std::min(newStepSize, getDeltaT()*2.0f); // For safety, limit how quickly dt can increase.
if (newStepSize > getDeltaT() && newStepSize < 1.2f*getDeltaT())
newStepSize = getDeltaT(); // Keeping dt constant between steps improves the behavior of the integrator.
if (newStepSize > maxStepSize)
newStepSize = maxStepSize;
RealOpenMM vstep = 0.5f*(newStepSize+getDeltaT()); // The time interval by which to advance the velocities
double vstep = 0.5f*(newStepSize+getDeltaT()); // The time interval by which to advance the velocities
setDeltaT(newStepSize);
for (int i = 0; i < numberOfAtoms; ++i) {
if (masses[i] != zero)
if (masses[i] != 0.0)
for (int j = 0; j < 3; ++j) {
RealOpenMM vPrime = velocities[i][j] + inverseMasses[i]*forces[i][j]*vstep;
double vPrime = velocities[i][j] + inverseMasses[i]*forces[i][j]*vstep;
xPrime[i][j] = atomCoordinates[i][j] + vPrime*getDeltaT();
}
}
......@@ -158,9 +141,9 @@ void ReferenceVariableVerletDynamics::update(const OpenMM::System& system, vecto
// Update the positions and velocities.
RealOpenMM velocityScale = one/getDeltaT();
double velocityScale = 1.0/getDeltaT();
for (int i = 0; i < numberOfAtoms; ++i) {
if (masses[i] != zero)
if (masses[i] != 0.0)
for (int j = 0; j < 3; ++j) {
velocities[i][j] = velocityScale*(xPrime[i][j] - atomCoordinates[i][j]);
atomCoordinates[i][j] = xPrime[i][j];
......
platforms/reference/src/SimTKReference/ReferenceVerletDynamics.cpp
View file @ a783b996
......@@ -45,18 +45,8 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */
ReferenceVerletDynamics::ReferenceVerletDynamics(int numberOfAtoms, RealOpenMM deltaT) :
ReferenceVerletDynamics::ReferenceVerletDynamics(int numberOfAtoms, double deltaT) :
ReferenceDynamics(numberOfAtoms, deltaT, 0.0) {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceVerletDynamics::ReferenceVerletDynamics";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
xPrime.resize(numberOfAtoms);
inverseMasses.resize(numberOfAtoms);
}
......@@ -68,13 +58,6 @@ ReferenceVerletDynamics::ReferenceVerletDynamics(int numberOfAtoms, RealOpenMM d
--------------------------------------------------------------------------------------- */
ReferenceVerletDynamics::~ReferenceVerletDynamics() {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceVerletDynamics::~ReferenceVerletDynamics";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
......@@ -90,19 +73,9 @@ ReferenceVerletDynamics::~ReferenceVerletDynamics() {
--------------------------------------------------------------------------------------- */
void ReferenceVerletDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
vector<RealVec>& velocities,
vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM tolerance) {
// ---------------------------------------------------------------------------------------
static const char* methodName = "\nReferenceVerletDynamics::update";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
void ReferenceVerletDynamics::update(const OpenMM::System& system, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities,
vector<Vec3>& forces, vector<double>& masses, double tolerance) {
// first-time-through initialization
int numberOfAtoms = system.getNumParticles();
......@@ -110,17 +83,17 @@ void ReferenceVerletDynamics::update(const OpenMM::System& system, vector<RealVe
// invert masses
for (int ii = 0; ii < numberOfAtoms; ii++) {
if (masses[ii] == zero)
inverseMasses[ii] = zero;
if (masses[ii] == 0.0)
inverseMasses[ii] = 0.0;
else
inverseMasses[ii] = one/masses[ii];
inverseMasses[ii] = 1.0/masses[ii];
}
}
// Perform the integration.
for (int i = 0; i < numberOfAtoms; ++i) {
if (masses[i] != zero)
if (masses[i] != 0.0)
for (int j = 0; j < 3; ++j) {
velocities[i][j] += inverseMasses[i]*forces[i][j]*getDeltaT();
xPrime[i][j] = atomCoordinates[i][j] + velocities[i][j]*getDeltaT();
......@@ -132,9 +105,9 @@ void ReferenceVerletDynamics::update(const OpenMM::System& system, vector<RealVe
// Update the positions and velocities.
RealOpenMM velocityScale = static_cast<RealOpenMM>(1.0/getDeltaT());
double velocityScale = static_cast<double>(1.0/getDeltaT());
for (int i = 0; i < numberOfAtoms; ++i) {
if (masses[i] != zero)
if (masses[i] != 0.0)
for (int j = 0; j < 3; ++j) {
velocities[i][j] = velocityScale*(xPrime[i][j] - atomCoordinates[i][j]);
atomCoordinates[i][j] = xPrime[i][j];
......
platforms/reference/src/SimTKReference/ReferenceVirtualSites.cpp
View file @ a783b996
......@@ -31,11 +31,12 @@
#include "ReferenceVirtualSites.h"
#include "openmm/VirtualSite.h"
#include <cmath>
using namespace OpenMM;
using namespace std;
void ReferenceVirtualSites::computePositions(const OpenMM::System& system, vector<OpenMM::RealVec>& atomCoordinates) {
void ReferenceVirtualSites::computePositions(const OpenMM::System& system, vector<OpenMM::Vec3>& atomCoordinates) {
for (int i = 0; i < system.getNumParticles(); i++)
if (system.isVirtualSite(i)) {
if (dynamic_cast<const TwoParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
......@@ -43,7 +44,7 @@ void ReferenceVirtualSites::computePositions(const OpenMM::System& system, vecto
const TwoParticleAverageSite& site = dynamic_cast<const TwoParticleAverageSite&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1);
RealOpenMM w1 = site.getWeight(0), w2 = site.getWeight(1);
double w1 = site.getWeight(0), w2 = site.getWeight(1);
atomCoordinates[i] = atomCoordinates[p1]*w1 + atomCoordinates[p2]*w2;
}
else if (dynamic_cast<const ThreeParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
......@@ -51,7 +52,7 @@ void ReferenceVirtualSites::computePositions(const OpenMM::System& system, vecto
const ThreeParticleAverageSite& site = dynamic_cast<const ThreeParticleAverageSite&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
RealOpenMM w1 = site.getWeight(0), w2 = site.getWeight(1), w3 = site.getWeight(2);
double w1 = site.getWeight(0), w2 = site.getWeight(1), w3 = site.getWeight(2);
atomCoordinates[i] = atomCoordinates[p1]*w1 + atomCoordinates[p2]*w2 + atomCoordinates[p3]*w3;
}
else if (dynamic_cast<const OutOfPlaneSite*>(&system.getVirtualSite(i)) != NULL) {
......@@ -59,10 +60,10 @@ void ReferenceVirtualSites::computePositions(const OpenMM::System& system, vecto
const OutOfPlaneSite& site = dynamic_cast<const OutOfPlaneSite&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
RealOpenMM w12 = site.getWeight12(), w13 = site.getWeight13(), wcross = site.getWeightCross();
RealVec v12 = atomCoordinates[p2]-atomCoordinates[p1];
RealVec v13 = atomCoordinates[p3]-atomCoordinates[p1];
RealVec cross = v12.cross(v13);
double w12 = site.getWeight12(), w13 = site.getWeight13(), wcross = site.getWeightCross();
Vec3 v12 = atomCoordinates[p2]-atomCoordinates[p1];
Vec3 v13 = atomCoordinates[p3]-atomCoordinates[p1];
Vec3 cross = v12.cross(v13);
atomCoordinates[i] = atomCoordinates[p1] + v12*w12 + v13*w13 + cross*wcross;
}
else if (dynamic_cast<const LocalCoordinatesSite*>(&system.getVirtualSite(i)) != NULL) {
......@@ -70,14 +71,14 @@ void ReferenceVirtualSites::computePositions(const OpenMM::System& system, vecto
const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
RealVec originWeights = site.getOriginWeights();
RealVec xWeights = site.getXWeights();
RealVec yWeights = site.getYWeights();
RealVec localPosition = site.getLocalPosition();
RealVec origin = atomCoordinates[p1]*originWeights[0] + atomCoordinates[p2]*originWeights[1] + atomCoordinates[p3]*originWeights[2];
RealVec xdir = atomCoordinates[p1]*xWeights[0] + atomCoordinates[p2]*xWeights[1] + atomCoordinates[p3]*xWeights[2];
RealVec ydir = atomCoordinates[p1]*yWeights[0] + atomCoordinates[p2]*yWeights[1] + atomCoordinates[p3]*yWeights[2];
RealVec zdir = xdir.cross(ydir);
Vec3 originWeights = site.getOriginWeights();
Vec3 xWeights = site.getXWeights();
Vec3 yWeights = site.getYWeights();
Vec3 localPosition = site.getLocalPosition();
Vec3 origin = atomCoordinates[p1]*originWeights[0] + atomCoordinates[p2]*originWeights[1] + atomCoordinates[p3]*originWeights[2];
Vec3 xdir = atomCoordinates[p1]*xWeights[0] + atomCoordinates[p2]*xWeights[1] + atomCoordinates[p3]*xWeights[2];
Vec3 ydir = atomCoordinates[p1]*yWeights[0] + atomCoordinates[p2]*yWeights[1] + atomCoordinates[p3]*yWeights[2];
Vec3 zdir = xdir.cross(ydir);
xdir /= sqrt(xdir.dot(xdir));
zdir /= sqrt(zdir.dot(zdir));
ydir = zdir.cross(xdir);
......@@ -87,16 +88,16 @@ void ReferenceVirtualSites::computePositions(const OpenMM::System& system, vecto
}
void ReferenceVirtualSites::distributeForces(const OpenMM::System& system, const vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& forces) {
void ReferenceVirtualSites::distributeForces(const OpenMM::System& system, const vector<OpenMM::Vec3>& atomCoordinates, vector<OpenMM::Vec3>& forces) {
for (int i = 0; i < system.getNumParticles(); i++)
if (system.isVirtualSite(i)) {
RealVec f = forces[i];
Vec3 f = forces[i];
if (dynamic_cast<const TwoParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
// A two particle average.
const TwoParticleAverageSite& site = dynamic_cast<const TwoParticleAverageSite&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1);
RealOpenMM w1 = site.getWeight(0), w2 = site.getWeight(1);
double w1 = site.getWeight(0), w2 = site.getWeight(1);
forces[p1] += f*w1;
forces[p2] += f*w2;
}
......@@ -105,7 +106,7 @@ void ReferenceVirtualSites::distributeForces(const OpenMM::System& system, const
const ThreeParticleAverageSite& site = dynamic_cast<const ThreeParticleAverageSite&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
RealOpenMM w1 = site.getWeight(0), w2 = site.getWeight(1), w3 = site.getWeight(2);
double w1 = site.getWeight(0), w2 = site.getWeight(1), w3 = site.getWeight(2);
forces[p1] += f*w1;
forces[p2] += f*w2;
forces[p3] += f*w3;
......@@ -115,15 +116,15 @@ void ReferenceVirtualSites::distributeForces(const OpenMM::System& system, const
const OutOfPlaneSite& site = dynamic_cast<const OutOfPlaneSite&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
RealOpenMM w12 = site.getWeight12(), w13 = site.getWeight13(), wcross = site.getWeightCross();
RealVec v12 = atomCoordinates[p2]-atomCoordinates[p1];
RealVec v13 = atomCoordinates[p3]-atomCoordinates[p1];
RealVec f2(w12*f[0] - wcross*v13[2]*f[1] + wcross*v13[1]*f[2],
wcross*v13[2]*f[0] + w12*f[1] - wcross*v13[0]*f[2],
-wcross*v13[1]*f[0] + wcross*v13[0]*f[1] + w12*f[2]);
RealVec f3(w13*f[0] + wcross*v12[2]*f[1] - wcross*v12[1]*f[2],
-wcross*v12[2]*f[0] + w13*f[1] + wcross*v12[0]*f[2],
wcross*v12[1]*f[0] - wcross*v12[0]*f[1] + w13*f[2]);
double w12 = site.getWeight12(), w13 = site.getWeight13(), wcross = site.getWeightCross();
Vec3 v12 = atomCoordinates[p2]-atomCoordinates[p1];
Vec3 v13 = atomCoordinates[p3]-atomCoordinates[p1];
Vec3 f2(w12*f[0] - wcross*v13[2]*f[1] + wcross*v13[1]*f[2],
wcross*v13[2]*f[0] + w12*f[1] - wcross*v13[0]*f[2],
-wcross*v13[1]*f[0] + wcross*v13[0]*f[1] + w12*f[2]);
Vec3 f3(w13*f[0] + wcross*v12[2]*f[1] - wcross*v12[1]*f[2],
-wcross*v12[2]*f[0] + w13*f[1] + wcross*v12[0]*f[2],
wcross*v12[1]*f[0] - wcross*v12[0]*f[1] + w13*f[2]);
forces[p1] += f-f2-f3;
forces[p2] += f2;
forces[p3] += f3;
......@@ -133,43 +134,43 @@ void ReferenceVirtualSites::distributeForces(const OpenMM::System& system, const
const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
RealVec originWeights = site.getOriginWeights();
RealVec wx = site.getXWeights();
RealVec wy = site.getYWeights();
RealVec localPosition = site.getLocalPosition();
RealVec xdir = atomCoordinates[p1]*wx[0] + atomCoordinates[p2]*wx[1] + atomCoordinates[p3]*wx[2];
RealVec ydir = atomCoordinates[p1]*wy[0] + atomCoordinates[p2]*wy[1] + atomCoordinates[p3]*wy[2];
RealVec zdir = xdir.cross(ydir);
RealOpenMM invNormXdir = 1.0/SQRT(xdir.dot(xdir));
RealOpenMM invNormZdir = 1.0/SQRT(zdir.dot(zdir));
RealVec dx = xdir*invNormXdir;
RealVec dz = zdir*invNormZdir;
RealVec dy = dz.cross(dx);
Vec3 originWeights = site.getOriginWeights();
Vec3 wx = site.getXWeights();
Vec3 wy = site.getYWeights();
Vec3 localPosition = site.getLocalPosition();
Vec3 xdir = atomCoordinates[p1]*wx[0] + atomCoordinates[p2]*wx[1] + atomCoordinates[p3]*wx[2];
Vec3 ydir = atomCoordinates[p1]*wy[0] + atomCoordinates[p2]*wy[1] + atomCoordinates[p3]*wy[2];
Vec3 zdir = xdir.cross(ydir);
double invNormXdir = 1.0/sqrt(xdir.dot(xdir));
double invNormZdir = 1.0/sqrt(zdir.dot(zdir));
Vec3 dx = xdir*invNormXdir;
Vec3 dz = zdir*invNormZdir;
Vec3 dy = dz.cross(dx);
// The derivatives for this case are very complicated. They were computed with SymPy then simplified by hand.
RealOpenMM t11 = (wx[0]*ydir[0]-wy[0]*xdir[0])*invNormZdir;
RealOpenMM t12 = (wx[0]*ydir[1]-wy[0]*xdir[1])*invNormZdir;
RealOpenMM t13 = (wx[0]*ydir[2]-wy[0]*xdir[2])*invNormZdir;
RealOpenMM t21 = (wx[1]*ydir[0]-wy[1]*xdir[0])*invNormZdir;
RealOpenMM t22 = (wx[1]*ydir[1]-wy[1]*xdir[1])*invNormZdir;
RealOpenMM t23 = (wx[1]*ydir[2]-wy[1]*xdir[2])*invNormZdir;
RealOpenMM t31 = (wx[2]*ydir[0]-wy[2]*xdir[0])*invNormZdir;
RealOpenMM t32 = (wx[2]*ydir[1]-wy[2]*xdir[1])*invNormZdir;
RealOpenMM t33 = (wx[2]*ydir[2]-wy[2]*xdir[2])*invNormZdir;
RealOpenMM sx1 = t13*dz[1]-t12*dz[2];
RealOpenMM sy1 = t11*dz[2]-t13*dz[0];
RealOpenMM sz1 = t12*dz[0]-t11*dz[1];
RealOpenMM sx2 = t23*dz[1]-t22*dz[2];
RealOpenMM sy2 = t21*dz[2]-t23*dz[0];
RealOpenMM sz2 = t22*dz[0]-t21*dz[1];
RealOpenMM sx3 = t33*dz[1]-t32*dz[2];
RealOpenMM sy3 = t31*dz[2]-t33*dz[0];
RealOpenMM sz3 = t32*dz[0]-t31*dz[1];
RealVec wxScaled = wx*invNormXdir;
RealVec fp1 = localPosition*f[0];
RealVec fp2 = localPosition*f[1];
RealVec fp3 = localPosition*f[2];
double t11 = (wx[0]*ydir[0]-wy[0]*xdir[0])*invNormZdir;
double t12 = (wx[0]*ydir[1]-wy[0]*xdir[1])*invNormZdir;
double t13 = (wx[0]*ydir[2]-wy[0]*xdir[2])*invNormZdir;
double t21 = (wx[1]*ydir[0]-wy[1]*xdir[0])*invNormZdir;
double t22 = (wx[1]*ydir[1]-wy[1]*xdir[1])*invNormZdir;
double t23 = (wx[1]*ydir[2]-wy[1]*xdir[2])*invNormZdir;
double t31 = (wx[2]*ydir[0]-wy[2]*xdir[0])*invNormZdir;
double t32 = (wx[2]*ydir[1]-wy[2]*xdir[1])*invNormZdir;
double t33 = (wx[2]*ydir[2]-wy[2]*xdir[2])*invNormZdir;
double sx1 = t13*dz[1]-t12*dz[2];
double sy1 = t11*dz[2]-t13*dz[0];
double sz1 = t12*dz[0]-t11*dz[1];
double sx2 = t23*dz[1]-t22*dz[2];
double sy2 = t21*dz[2]-t23*dz[0];
double sz2 = t22*dz[0]-t21*dz[1];
double sx3 = t33*dz[1]-t32*dz[2];
double sy3 = t31*dz[2]-t33*dz[0];
double sz3 = t32*dz[0]-t31*dz[1];
Vec3 wxScaled = wx*invNormXdir;
Vec3 fp1 = localPosition*f[0];
Vec3 fp2 = localPosition*f[1];
Vec3 fp3 = localPosition*f[2];
forces[p1][0] += fp1[0]*wxScaled[0]*(1-dx[0]*dx[0]) + fp1[2]*(dz[0]*sx1 ) + fp1[1]*((-dx[0]*dy[0] )*wxScaled[0] + dy[0]*sx1 - dx[1]*t12 - dx[2]*t13) + f[0]*originWeights[0];
forces[p1][1] += fp1[0]*wxScaled[0]*( -dx[0]*dx[1]) + fp1[2]*(dz[0]*sy1+t13) + fp1[1]*((-dx[1]*dy[0]-dz[2])*wxScaled[0] + dy[0]*sy1 + dx[1]*t11);
forces[p1][2] += fp1[0]*wxScaled[0]*( -dx[0]*dx[2]) + fp1[2]*(dz[0]*sz1-t12) + fp1[1]*((-dx[2]*dy[0]+dz[1])*wxScaled[0] + dy[0]*sz1 + dx[2]*t11);
......
platforms/reference/src/SimTKReference/SimTKOpenMMUtilities.cpp
View file @ a783b996
......@@ -40,12 +40,12 @@ using namespace OpenMM;
uint32_t SimTKOpenMMUtilities::_randomNumberSeed = 0;
bool SimTKOpenMMUtilities::_randomInitialized = false;
bool SimTKOpenMMUtilities::nextGaussianIsValid = false;
RealOpenMM SimTKOpenMMUtilities::nextGaussian = 0;
double SimTKOpenMMUtilities::nextGaussian = 0;
OpenMM_SFMT::SFMT SimTKOpenMMUtilities::sfmt;
/* ---------------------------------------------------------------------------------------
Allocate 1D RealOpenMM array (Simbios)
Allocate 1D double array (Simbios)
array[i]
......@@ -59,27 +59,18 @@ OpenMM_SFMT::SFMT SimTKOpenMMUtilities::sfmt;
--------------------------------------------------------------------------------------- */
RealOpenMM* SimTKOpenMMUtilities::allocateOneDRealOpenMMArray(int iSize, RealOpenMM* array1D,
int initialize, RealOpenMM initialValue,
double* SimTKOpenMMUtilities::allocateOneDRealOpenMMArray(int iSize, double* array1D,
int initialize, double initialValue,
const std::string& idString) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nSimTKOpenMMUtilities::allocate1DRealOpenMMArray";
static const RealOpenMM zero = 0.0;
// ---------------------------------------------------------------------------------------
if (array1D == NULL) {
array1D = new RealOpenMM[iSize];
array1D = new double[iSize];
}
if (initialize) {
if (initialValue == zero) {
memset(array1D, 0, iSize*sizeof(RealOpenMM));
if (initialValue == 0.0) {
memset(array1D, 0, iSize*sizeof(double));
} else {
for (int ii = 0; ii < iSize; ii++) {
array1D[ii] = initialValue;
......@@ -92,7 +83,7 @@ RealOpenMM* SimTKOpenMMUtilities::allocateOneDRealOpenMMArray(int iSize, RealOpe
/* ---------------------------------------------------------------------------------------
Allocate 2D RealOpenMM array (Simbios)
Allocate 2D double array (Simbios)
array[i][j]
......@@ -107,23 +98,16 @@ RealOpenMM* SimTKOpenMMUtilities::allocateOneDRealOpenMMArray(int iSize, RealOpe
--------------------------------------------------------------------------------------- */
RealOpenMM** SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray(int iSize, int jSize, RealOpenMM** array2D,
int initialize, RealOpenMM initialValue,
double** SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray(int iSize, int jSize, double** array2D,
int initialize, double initialValue,
const std::string& idString) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nSimTKOpenMMUtilities::allocate2DRealOpenMMArray";
// ---------------------------------------------------------------------------------------
if (array2D == NULL) {
array2D = new RealOpenMM*[iSize];
array2D = new double*[iSize];
std::string blockString = idString;
blockString.append("Block");
RealOpenMM* block = new RealOpenMM[jSize*iSize];
double* block = new double[jSize*iSize];
for (int ii = 0; ii < iSize; ii++) {
array2D[ii] = block;
......@@ -140,7 +124,7 @@ RealOpenMM** SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray(int iSize, int jS
/* ---------------------------------------------------------------------------------------
Free 2D RealOpenMM array (Simbios)
Free 2D double array (Simbios)
array[i][j]
......@@ -149,14 +133,7 @@ RealOpenMM** SimTKOpenMMUtilities::allocateTwoDRealOpenMMArray(int iSize, int jS
--------------------------------------------------------------------------------------- */
void SimTKOpenMMUtilities::freeTwoDRealOpenMMArray(RealOpenMM** array2D, const std::string& idString) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nSimTKOpenMMUtilities::freeTwoDRealOpenMMArray";
// ---------------------------------------------------------------------------------------
void SimTKOpenMMUtilities::freeTwoDRealOpenMMArray(double** array2D, const std::string& idString) {
if (array2D != NULL) {
std::string blockString = idString;
......@@ -169,7 +146,7 @@ void SimTKOpenMMUtilities::freeTwoDRealOpenMMArray(RealOpenMM** array2D, const s
/* ---------------------------------------------------------------------------------------
Free 1D RealOpenMM array (Simbios)
Free 1D double array (Simbios)
array[i]
......@@ -178,14 +155,7 @@ void SimTKOpenMMUtilities::freeTwoDRealOpenMMArray(RealOpenMM** array2D, const s
--------------------------------------------------------------------------------------- */
void SimTKOpenMMUtilities::freeOneDRealOpenMMArray(RealOpenMM* array1D, const std::string& idString) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nSimTKOpenMMUtilities::freeOneDRealOpenMMArray";
// ---------------------------------------------------------------------------------------
void SimTKOpenMMUtilities::freeOneDRealOpenMMArray(double* array1D, const std::string& idString) {
if (array1D != NULL) {
delete[] array1D;
}
......@@ -193,7 +163,7 @@ void SimTKOpenMMUtilities::freeOneDRealOpenMMArray(RealOpenMM* array1D, const st
/* ---------------------------------------------------------------------------------------
Initialize 2D RealOpenMM array (Simbios)
Initialize 2D double array (Simbios)
array[i][j]
......@@ -205,15 +175,8 @@ void SimTKOpenMMUtilities::freeOneDRealOpenMMArray(RealOpenMM* array1D, const st
--------------------------------------------------------------------------------------- */
void SimTKOpenMMUtilities::initialize2DRealOpenMMArray(int iSize, int jSize,
RealOpenMM** array2D,
RealOpenMM initialValue) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nSimTKOpenMMUtilities::initialize2DRealOpenMMArray";
// ---------------------------------------------------------------------------------------
double** array2D,
double initialValue) {
bool useMemset;
bool useMemsetSingleBlock;
......@@ -231,10 +194,10 @@ void SimTKOpenMMUtilities::initialize2DRealOpenMMArray(int iSize, int jSize,
if (useMemset) {
if (useMemsetSingleBlock) {
memset(array2D[0], 0, iSize*jSize*sizeof(RealOpenMM));
memset(array2D[0], 0, iSize*jSize*sizeof(double));
} else {
for (int ii = 0; ii < iSize; ii++) {
memset(array2D[ii], 0, jSize*sizeof(RealOpenMM));
memset(array2D[ii], 0, jSize*sizeof(double));
}
}
} else {
......@@ -260,16 +223,9 @@ void SimTKOpenMMUtilities::initialize2DRealOpenMMArray(int iSize, int jSize,
--------------------------------------------------------------------------------------- */
void SimTKOpenMMUtilities::crossProductVector3(RealOpenMM* vectorX,
RealOpenMM* vectorY,
RealOpenMM* vectorZ) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nSimTKOpenMMUtilities::crossProductVector3";
// ---------------------------------------------------------------------------------------
void SimTKOpenMMUtilities::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];
vectorZ[2] = vectorX[0]*vectorY[1] - vectorX[1]*vectorY[0];
......@@ -285,7 +241,7 @@ void SimTKOpenMMUtilities::crossProductVector3(RealOpenMM* vectorX,
--------------------------------------------------------------------------------------- */
RealOpenMM SimTKOpenMMUtilities::getNormallyDistributedRandomNumber() {
double SimTKOpenMMUtilities::getNormallyDistributedRandomNumber() {
if (nextGaussianIsValid) {
nextGaussianIsValid = false;
return nextGaussian;
......@@ -298,13 +254,13 @@ RealOpenMM SimTKOpenMMUtilities::getNormallyDistributedRandomNumber() {
// Use the polar form of the Box-Muller transformation to generate two Gaussian random numbers.
RealOpenMM x, y, r2;
double x, y, r2;
do {
x = static_cast<RealOpenMM>(2.0 * genrand_real2(sfmt) - 1.0);
y = static_cast<RealOpenMM>(2.0 * genrand_real2(sfmt) - 1.0);
x = 2.0 * genrand_real2(sfmt) - 1.0;
y = 2.0 * genrand_real2(sfmt) - 1.0;
r2 = x*x + y*y;
} while (r2 >= 1.0 || r2 == 0.0);
RealOpenMM multiplier = static_cast<RealOpenMM>(sqrt((-2.0*log(r2))/r2));
double multiplier = sqrt((-2.0*log(r2))/r2);
nextGaussian = y*multiplier;
nextGaussianIsValid = true;
return x*multiplier;
......@@ -318,13 +274,13 @@ RealOpenMM SimTKOpenMMUtilities::getNormallyDistributedRandomNumber() {
--------------------------------------------------------------------------------------- */
RealOpenMM SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() {
double SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() {
if (!_randomInitialized) {
init_gen_rand(_randomNumberSeed, sfmt);
_randomInitialized = true;
nextGaussianIsValid = false;
}
RealOpenMM value = static_cast<RealOpenMM>(genrand_real2(sfmt));
double value = genrand_real2(sfmt);
return value;
}
......@@ -337,13 +293,6 @@ RealOpenMM SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() {
--------------------------------------------------------------------------------------- */
uint32_t SimTKOpenMMUtilities::getRandomNumberSeed() {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceDynamics::getRandomNumberSeed";
// ---------------------------------------------------------------------------------------
return _randomNumberSeed;
}
......@@ -356,13 +305,6 @@ uint32_t SimTKOpenMMUtilities::getRandomNumberSeed() {
--------------------------------------------------------------------------------------- */
void SimTKOpenMMUtilities::setRandomNumberSeed(uint32_t seed) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceDynamics::setRandomNumberSeed";
// ---------------------------------------------------------------------------------------
// If the seed is 0, use a unique seed
if (seed == 0)
_randomNumberSeed = (uint32_t) osrngseed();
......@@ -376,7 +318,7 @@ void SimTKOpenMMUtilities::createCheckpoint(std::ostream& stream) {
stream.write((char*) &_randomInitialized, sizeof(bool));
if (_randomInitialized) {
stream.write((char*) &nextGaussianIsValid, sizeof(bool));
stream.write((char*) &nextGaussian, sizeof(RealOpenMM));
stream.write((char*) &nextGaussian, sizeof(double));
sfmt.createCheckpoint(stream);
}
}
......@@ -389,7 +331,7 @@ void SimTKOpenMMUtilities::loadCheckpoint(std::istream& stream) {
if (!prevInitialized)
init_gen_rand(0, sfmt);
stream.read((char*) &nextGaussianIsValid, sizeof(bool));
stream.read((char*) &nextGaussian, sizeof(RealOpenMM));
stream.read((char*) &nextGaussian, sizeof(double));
sfmt.loadCheckpoint(stream);
}
}
platforms/reference/src/SimTKReference/fftpack.cpp
View file @ a783b996
......@@ -37,7 +37,7 @@ struct fftpack
int n; /**< Number of points in this dimension. */
int ifac[15]; /**< 15 bytes needed for cfft and rfft */
struct fftpack * next; /**< Pointer to next dimension, or NULL. */
RealOpenMM * work; /**< 1st 4n reserved for cfft, 1st 2n for rfft */
double * work; /**< 1st 4n reserved for cfft, 1st 2n for rfft */
};
......@@ -48,15 +48,15 @@ struct fftpack
static void
fftpack_passf2(int ido,
int l1,
RealOpenMM cc[],
RealOpenMM ch[],
RealOpenMM wa1[],
int isign)
fftpack_passf2(int ido,
int l1,
double cc[],
double ch[],
double wa1[],
int isign)
{
int i, k, ah, ac;
RealOpenMM ti2, tr2;
double ti2, tr2;
if (ido <= 2)
{
......@@ -92,19 +92,19 @@ fftpack_passf2(int ido,
static void
fftpack_passf3(int ido,
int l1,
RealOpenMM cc[],
RealOpenMM ch[],
RealOpenMM wa1[],
RealOpenMM wa2[],
int isign)
fftpack_passf3(int ido,
int l1,
double cc[],
double ch[],
double wa1[],
double wa2[],
int isign)
{
const RealOpenMM taur = (RealOpenMM)-0.5;
const RealOpenMM taui = (RealOpenMM)0.866025403784439;
const double taur = -0.5;
const double taui = 0.866025403784439;
int i, k, ac, ah;
RealOpenMM ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
if (ido == 2)
{
......@@ -159,17 +159,17 @@ fftpack_passf3(int ido,
static void
fftpack_passf4(int ido,
int l1,
RealOpenMM cc[],
RealOpenMM ch[],
RealOpenMM wa1[],
RealOpenMM wa2[],
RealOpenMM wa3[],
int isign)
fftpack_passf4(int ido,
int l1,
double cc[],
double ch[],
double wa1[],
double wa2[],
double wa3[],
int isign)
{
int i, k, ac, ah;
RealOpenMM ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
if (ido == 2)
{
......@@ -232,23 +232,23 @@ fftpack_passf4(int ido,
static void
fftpack_passf5(int ido,
int l1,
RealOpenMM cc[],
RealOpenMM ch[],
RealOpenMM wa1[],
RealOpenMM wa2[],
RealOpenMM wa3[],
RealOpenMM wa4[],
int isign)
fftpack_passf5(int ido,
int l1,
double cc[],
double ch[],
double wa1[],
double wa2[],
double wa3[],
double wa4[],
int isign)
{
const RealOpenMM tr11 = (RealOpenMM)0.309016994374947;
const RealOpenMM ti11 = (RealOpenMM)0.951056516295154;
const RealOpenMM tr12 = (RealOpenMM)-0.809016994374947;
const RealOpenMM ti12 = (RealOpenMM)0.587785252292473;
const double tr11 = 0.309016994374947;
const double ti11 = 0.951056516295154;
const double tr12 = -0.809016994374947;
const double ti12 = 0.587785252292473;
int i, k, ac, ah;
RealOpenMM ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3,
ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
if (ido == 2)
......@@ -334,18 +334,18 @@ fftpack_passf5(int ido,
static void
fftpack_passf(int * nac,
int ido,
int ip,
int l1,
int idl1,
RealOpenMM cc[],
RealOpenMM ch[],
RealOpenMM wa[],
int isign)
fftpack_passf(int* nac,
int ido,
int ip,
int l1,
int idl1,
double cc[],
double ch[],
double wa[],
int isign)
{
int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, nt, idj, idl, inc,idp;
RealOpenMM wai, war;
double wai, war;
idot = ido / 2;
nt = ip*idl1;
......@@ -495,17 +495,17 @@ fftpack_passf(int * nac,
static void
fftpack_cfftf1(int n,
RealOpenMM c[],
RealOpenMM ch[],
RealOpenMM wa[],
int ifac[15],
int isign)
fftpack_cfftf1(int n,
double c[],
double ch[],
double wa[],
int ifac[15],
int isign)
{
int idot, i;
int k1, l1, l2;
int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1;
RealOpenMM *cinput, *coutput;
double *cinput, *coutput;
nf = ifac[1];
na = 0;
l1 = 1;
......@@ -606,12 +606,12 @@ startloop:
static void
fftpack_cffti1(int n,
RealOpenMM wa[],
int ifac[15])
fftpack_cffti1(int n,
double wa[],
int ifac[15])
{
const RealOpenMM twopi = (RealOpenMM)6.28318530717959;
RealOpenMM arg, argh, argld, fi;
const double twopi = 6.28318530717959;
double arg, argh, argld, fi;
int idot, i, j;
int i1, k1, l1, l2;
int ld, ii, nf, ip;
......@@ -619,7 +619,7 @@ fftpack_cffti1(int n,
fftpack_factorize(n,ifac);
nf = ifac[1];
argh = twopi/(RealOpenMM)n;
argh = twopi/n;
i = 1;
l1 = 1;
for (k1=1; k1<=nf; k1++)
......@@ -783,7 +783,7 @@ fftpack_init_1d(fftpack_t * pfft,
fft->n = nx;
/* Need 4*n storage for 1D complex FFT */
if ((fft->work = (RealOpenMM *)malloc(sizeof(RealOpenMM)*(4*nx))) == NULL)
if ((fft->work = (double *)malloc(sizeof(double)*(4*nx))) == NULL)
{
free(fft);
return ENOMEM;
......@@ -860,7 +860,7 @@ fftpack_init_3d(fftpack_t * pfft,
/* Need 4*nx storage for 1D complex FFT.
*/
if ((fft->work = (RealOpenMM *)malloc(sizeof(RealOpenMM)*(4*nx))) == NULL)
if ((fft->work = (double *)malloc(sizeof(double)*(4*nx))) == NULL)
{
free(fft);
return ENOMEM;
......@@ -887,16 +887,16 @@ fftpack_exec_1d (fftpack_t fft,
t_complex * in_data,
t_complex * out_data)
{
int i,n;
RealOpenMM * p1;
RealOpenMM * p2;
int i, n;
double* p1;
double* p2;
n=fft->n;
if (n==1)
{
p1 = (RealOpenMM *)in_data;
p2 = (RealOpenMM *)out_data;
p1 = (double *)in_data;
p2 = (double *)out_data;
p2[0] = p1[0];
p2[1] = p1[1];
}
......@@ -906,10 +906,10 @@ fftpack_exec_1d (fftpack_t fft,
*/
if (in_data != out_data)
{
p1 = (RealOpenMM *)in_data;
p2 = (RealOpenMM *)out_data;
p1 = (double *)in_data;
p2 = (double *)out_data;
/* n complex = 2*n RealOpenMM elements */
/* n complex = 2*n double elements */
for (i=0;i<2*n;i++)
{
p2[i] = p1[i];
......@@ -922,11 +922,11 @@ fftpack_exec_1d (fftpack_t fft,
if (dir == FFTPACK_FORWARD)
{
fftpack_cfftf1(n,(RealOpenMM *)out_data,fft->work+2*n,fft->work,fft->ifac, -1);
fftpack_cfftf1(n,(double *)out_data,fft->work+2*n,fft->work,fft->ifac, -1);
}
else if (dir == FFTPACK_BACKWARD)
{
fftpack_cfftf1(n,(RealOpenMM *)out_data,fft->work+2*n,fft->work,fft->ifac, 1);
fftpack_cfftf1(n,(double *)out_data,fft->work+2*n,fft->work,fft->ifac, 1);
}
else
{
......
platforms/reference/tests/TestReferenceNeighborList.cpp
View file @ a783b996
......@@ -41,14 +41,14 @@ using namespace OpenMM;
void testNeighborList()
{
vector<RealVec> particleList(2);
particleList[0] = RealVec(13.6, 0, 0);
particleList[1] = RealVec(0, 0, 0);
vector<Vec3> particleList(2);
particleList[0] = Vec3(13.6, 0, 0);
particleList[1] = Vec3(0, 0, 0);
vector<set<int> > exclusions(2);
NeighborList neighborList;
RealVec boxVectors[3];
Vec3 boxVectors[3];
computeNeighborListNaive(neighborList, 2, particleList, exclusions, boxVectors, false, 13.7, 0.01);
assert(neighborList.size() == 1);
......@@ -62,15 +62,15 @@ void testNeighborList()
assert(neighborList.size() == 0);
}
double distance2(RealVec& pos1, RealVec& pos2, const RealVec* periodicBoxVectors) {
RealVec diff = pos1-pos2;
double distance2(Vec3& pos1, Vec3& pos2, const Vec3* periodicBoxVectors) {
Vec3 diff = pos1-pos2;
diff -= periodicBoxVectors[2]*floor(diff[2]/periodicBoxVectors[2][2]+0.5);
diff -= periodicBoxVectors[1]*floor(diff[1]/periodicBoxVectors[1][1]+0.5);
diff -= periodicBoxVectors[0]*floor(diff[0]/periodicBoxVectors[0][0]+0.5);
return diff.dot(diff);
}
void verifyNeighborList(NeighborList& list, int numParticles, vector<RealVec>& positions, const RealVec* periodicBoxVectors, double cutoff) {
void verifyNeighborList(NeighborList& list, int numParticles, vector<Vec3>& positions, const Vec3* periodicBoxVectors, double cutoff) {
for (int i = 0; i < (int) list.size(); i++) {
int particle1 = list[i].first;
int particle2 = list[i].second;
......@@ -87,18 +87,18 @@ void verifyNeighborList(NeighborList& list, int numParticles, vector<RealVec>& p
void testPeriodic() {
const int numParticles = 100;
const double cutoff = 3.0;
RealVec periodicBoxVectors[3];
periodicBoxVectors[0] = RealVec(20, 0, 0);
periodicBoxVectors[1] = RealVec(0, 15, 0);
periodicBoxVectors[2] = RealVec(0, 0, 22);
vector<RealVec> particleList(numParticles);
Vec3 periodicBoxVectors[3];
periodicBoxVectors[0] = Vec3(20, 0, 0);
periodicBoxVectors[1] = Vec3(0, 15, 0);
periodicBoxVectors[2] = Vec3(0, 0, 22);
vector<Vec3> particleList(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i <numParticles; i++) {
particleList[i][0] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[0][0]*3);
particleList[i][1] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[1][1]*3);
particleList[i][2] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[2][2]*3);
particleList[i][0] = genrand_real2(sfmt)*periodicBoxVectors[0][0]*3;
particleList[i][1] = genrand_real2(sfmt)*periodicBoxVectors[1][1]*3;
particleList[i][2] = genrand_real2(sfmt)*periodicBoxVectors[2][2]*3;
}
vector<set<int> > exclusions(numParticles);
NeighborList neighborList;
......@@ -111,18 +111,18 @@ void testPeriodic() {
void testTriclinic() {
const int numParticles = 1000;
const double cutoff = 3.0;
RealVec periodicBoxVectors[3];
periodicBoxVectors[0] = RealVec(20, 0, 0);
periodicBoxVectors[1] = RealVec(5, 15, 0);
periodicBoxVectors[2] = RealVec(-3, -7, 22);
vector<RealVec> particleList(numParticles);
Vec3 periodicBoxVectors[3];
periodicBoxVectors[0] = Vec3(20, 0, 0);
periodicBoxVectors[1] = Vec3(5, 15, 0);
periodicBoxVectors[2] = Vec3(-3, -7, 22);
vector<Vec3> particleList(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i <numParticles; i++) {
particleList[i][0] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[0][0]*3);
particleList[i][1] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[1][1]*3);
particleList[i][2] = (RealOpenMM) (genrand_real2(sfmt)*periodicBoxVectors[2][2]*3);
particleList[i][0] = genrand_real2(sfmt)*periodicBoxVectors[0][0]*3;
particleList[i][1] = genrand_real2(sfmt)*periodicBoxVectors[1][1]*3;
particleList[i][2] = genrand_real2(sfmt)*periodicBoxVectors[2][2]*3;
}
vector<set<int> > exclusions(numParticles);
NeighborList neighborList;
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ a783b996
......@@ -54,29 +54,29 @@
using namespace OpenMM;
using namespace std;
static vector<RealVec>& extractPositions(ContextImpl& context) {
static vector<Vec3>& extractPositions(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<RealVec>*) data->positions);
return *((vector<Vec3>*) data->positions);
}
static vector<RealVec>& extractVelocities(ContextImpl& context) {
static vector<Vec3>& extractVelocities(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<RealVec>*) data->velocities);
return *((vector<Vec3>*) data->velocities);
}
static vector<RealVec>& extractForces(ContextImpl& context) {
static vector<Vec3>& extractForces(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<RealVec>*) data->forces);
return *((vector<Vec3>*) data->forces);
}
static RealVec& extractBoxSize(ContextImpl& context) {
static Vec3& extractBoxSize(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *(RealVec*) data->periodicBoxSize;
return *(Vec3*) data->periodicBoxSize;
}
static RealVec* extractBoxVectors(ContextImpl& context) {
static Vec3* extractBoxVectors(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return (RealVec*) data->periodicBoxVectors;
return (Vec3*) data->periodicBoxVectors;
}
// ***************************************************************************
......@@ -99,23 +99,23 @@ void ReferenceCalcAmoebaBondForceKernel::initialize(const System& system, const
particle1.push_back(particle1Index);
particle2.push_back(particle2Index);
length.push_back(static_cast<RealOpenMM>(lengthValue));
kQuadratic.push_back(static_cast<RealOpenMM>(kValue));
length.push_back(static_cast<double>(lengthValue));
kQuadratic.push_back(kValue);
}
globalBondCubic = static_cast<RealOpenMM>(force.getAmoebaGlobalBondCubic());
globalBondQuartic = static_cast<RealOpenMM>(force.getAmoebaGlobalBondQuartic());
globalBondCubic = force.getAmoebaGlobalBondCubic();
globalBondQuartic = force.getAmoebaGlobalBondQuartic();
usePeriodic = force.usesPeriodicBoundaryConditions();
}
double ReferenceCalcAmoebaBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
AmoebaReferenceBondForce amoebaReferenceBondForce;
if (usePeriodic)
amoebaReferenceBondForce.setPeriodic(extractBoxVectors(context));
RealOpenMM energy = amoebaReferenceBondForce.calculateForceAndEnergy(numBonds, posData, particle1, particle2, length, kQuadratic,
globalBondCubic, globalBondQuartic,
forceData);
double energy = amoebaReferenceBondForce.calculateForceAndEnergy(numBonds, posData, particle1, particle2, length, kQuadratic,
globalBondCubic, globalBondQuartic,
forceData);
return static_cast<double>(energy);
}
......@@ -131,8 +131,8 @@ void ReferenceCalcAmoebaBondForceKernel::copyParametersToContext(ContextImpl& co
force.getBondParameters(i, particle1Index, particle2Index, lengthValue, kValue);
if (particle1Index != particle1[i] || particle2Index != particle2[i])
throw OpenMMException("updateParametersInContext: The set of particles in a bond has changed");
length[i] = (RealOpenMM) lengthValue;
kQuadratic[i] = (RealOpenMM) kValue;
length[i] = lengthValue;
kQuadratic[i] = kValue;
}
}
......@@ -156,23 +156,23 @@ void ReferenceCalcAmoebaAngleForceKernel::initialize(const System& system, const
particle1.push_back(particle1Index);
particle2.push_back(particle2Index);
particle3.push_back(particle3Index);
angle.push_back(static_cast<RealOpenMM>(angleValue));
kQuadratic.push_back(static_cast<RealOpenMM>(k));
angle.push_back(angleValue);
kQuadratic.push_back(k);
}
globalAngleCubic = static_cast<RealOpenMM>(force.getAmoebaGlobalAngleCubic());
globalAngleQuartic = static_cast<RealOpenMM>(force.getAmoebaGlobalAngleQuartic());
globalAnglePentic = static_cast<RealOpenMM>(force.getAmoebaGlobalAnglePentic());
globalAngleSextic = static_cast<RealOpenMM>(force.getAmoebaGlobalAngleSextic());
globalAngleCubic = force.getAmoebaGlobalAngleCubic();
globalAngleQuartic = force.getAmoebaGlobalAngleQuartic();
globalAnglePentic = force.getAmoebaGlobalAnglePentic();
globalAngleSextic = force.getAmoebaGlobalAngleSextic();
usePeriodic = force.usesPeriodicBoundaryConditions();
}
double ReferenceCalcAmoebaAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
AmoebaReferenceAngleForce amoebaReferenceAngleForce;
if (usePeriodic)
amoebaReferenceAngleForce.setPeriodic(extractBoxVectors(context));
RealOpenMM energy = amoebaReferenceAngleForce.calculateForceAndEnergy(numAngles,
double energy = amoebaReferenceAngleForce.calculateForceAndEnergy(numAngles,
posData, particle1, particle2, particle3, angle, kQuadratic, globalAngleCubic, globalAngleQuartic, globalAnglePentic, globalAngleSextic, forceData);
return static_cast<double>(energy);
}
......@@ -189,8 +189,8 @@ void ReferenceCalcAmoebaAngleForceKernel::copyParametersToContext(ContextImpl& c
force.getAngleParameters(i, particle1Index, particle2Index, particle3Index, angleValue, k);
if (particle1Index != particle1[i] || particle2Index != particle2[i] || particle3Index != particle3[i])
throw OpenMMException("updateParametersInContext: The set of particles in an angle has changed");
angle[i] = (RealOpenMM) angleValue;
kQuadratic[i] = (RealOpenMM) k;
angle[i] = angleValue;
kQuadratic[i] = k;
}
}
......@@ -212,26 +212,26 @@ void ReferenceCalcAmoebaInPlaneAngleForceKernel::initialize(const System& system
particle2.push_back(particle2Index);
particle3.push_back(particle3Index);
particle4.push_back(particle4Index);
angle.push_back(static_cast<RealOpenMM>(angleValue));
kQuadratic.push_back(static_cast<RealOpenMM>(k));
angle.push_back(angleValue);
kQuadratic.push_back(k);
}
globalInPlaneAngleCubic = static_cast<RealOpenMM>(force.getAmoebaGlobalInPlaneAngleCubic());
globalInPlaneAngleQuartic = static_cast<RealOpenMM>(force.getAmoebaGlobalInPlaneAngleQuartic());
globalInPlaneAnglePentic = static_cast<RealOpenMM>(force.getAmoebaGlobalInPlaneAnglePentic());
globalInPlaneAngleSextic = static_cast<RealOpenMM>(force.getAmoebaGlobalInPlaneAngleSextic());
globalInPlaneAngleCubic = force.getAmoebaGlobalInPlaneAngleCubic();
globalInPlaneAngleQuartic = force.getAmoebaGlobalInPlaneAngleQuartic();
globalInPlaneAnglePentic = force.getAmoebaGlobalInPlaneAnglePentic();
globalInPlaneAngleSextic = force.getAmoebaGlobalInPlaneAngleSextic();
usePeriodic = force.usesPeriodicBoundaryConditions();
}
double ReferenceCalcAmoebaInPlaneAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
AmoebaReferenceInPlaneAngleForce amoebaReferenceInPlaneAngleForce;
if (usePeriodic)
amoebaReferenceInPlaneAngleForce.setPeriodic(extractBoxVectors(context));
RealOpenMM energy = amoebaReferenceInPlaneAngleForce.calculateForceAndEnergy(numAngles, posData, particle1, particle2, particle3, particle4,
angle, kQuadratic, globalInPlaneAngleCubic, globalInPlaneAngleQuartic,
globalInPlaneAnglePentic, globalInPlaneAngleSextic, forceData);
double energy = amoebaReferenceInPlaneAngleForce.calculateForceAndEnergy(numAngles, posData, particle1, particle2, particle3, particle4,
angle, kQuadratic, globalInPlaneAngleCubic, globalInPlaneAngleQuartic,
globalInPlaneAnglePentic, globalInPlaneAngleSextic, forceData);
return static_cast<double>(energy);
}
......@@ -247,8 +247,8 @@ void ReferenceCalcAmoebaInPlaneAngleForceKernel::copyParametersToContext(Context
force.getAngleParameters(i, particle1Index, particle2Index, particle3Index, particle4Index, angleValue, k);
if (particle1Index != particle1[i] || particle2Index != particle2[i] || particle3Index != particle3[i] || particle4Index != particle4[i])
throw OpenMMException("updateParametersInContext: The set of particles in an angle has changed");
angle[i] = (RealOpenMM) angleValue;
kQuadratic[i] = (RealOpenMM) k;
angle[i] = angleValue;
kQuadratic[i] = k;
}
}
......@@ -273,18 +273,18 @@ void ReferenceCalcAmoebaPiTorsionForceKernel::initialize(const System& system, c
particle4.push_back(particle4Index);
particle5.push_back(particle5Index);
particle6.push_back(particle6Index);
kTorsion.push_back(static_cast<RealOpenMM>(kTorsionParameter));
kTorsion.push_back(kTorsionParameter);
}
usePeriodic = force.usesPeriodicBoundaryConditions();
}
double ReferenceCalcAmoebaPiTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
AmoebaReferencePiTorsionForce amoebaReferencePiTorsionForce;
if (usePeriodic)
amoebaReferencePiTorsionForce.setPeriodic(extractBoxVectors(context));
RealOpenMM energy = amoebaReferencePiTorsionForce.calculateForceAndEnergy(numPiTorsions, posData, particle1, particle2,
double energy = amoebaReferencePiTorsionForce.calculateForceAndEnergy(numPiTorsions, posData, particle1, particle2,
particle3, particle4, particle5, particle6,
kTorsion, forceData);
return static_cast<double>(energy);
......@@ -303,7 +303,7 @@ void ReferenceCalcAmoebaPiTorsionForceKernel::copyParametersToContext(ContextImp
if (particle1Index != particle1[i] || particle2Index != particle2[i] || particle3Index != particle3[i] ||
particle4Index != particle4[i] || particle5Index != particle5[i] || particle6Index != particle6[i])
throw OpenMMException("updateParametersInContext: The set of particles in a torsion has changed");
kTorsion[i] = (RealOpenMM) kTorsionParameter;
kTorsion[i] = kTorsionParameter;
}
}
......@@ -324,22 +324,22 @@ void ReferenceCalcAmoebaStretchBendForceKernel::initialize(const System& system,
particle1.push_back(particle1Index);
particle2.push_back(particle2Index);
particle3.push_back(particle3Index);
lengthABParameters.push_back(static_cast<RealOpenMM>(lengthAB));
lengthCBParameters.push_back(static_cast<RealOpenMM>(lengthCB));
angleParameters.push_back(static_cast<RealOpenMM>(angle));
k1Parameters.push_back(static_cast<RealOpenMM>(k1));
k2Parameters.push_back(static_cast<RealOpenMM>(k2));
lengthABParameters.push_back(lengthAB);
lengthCBParameters.push_back(lengthCB);
angleParameters.push_back(angle);
k1Parameters.push_back(k1);
k2Parameters.push_back(k2);
}
usePeriodic = force.usesPeriodicBoundaryConditions();
}
double ReferenceCalcAmoebaStretchBendForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
AmoebaReferenceStretchBendForce amoebaReferenceStretchBendForce;
if (usePeriodic)
amoebaReferenceStretchBendForce.setPeriodic(extractBoxVectors(context));
RealOpenMM energy = amoebaReferenceStretchBendForce.calculateForceAndEnergy(numStretchBends, posData, particle1, particle2, particle3,
double energy = amoebaReferenceStretchBendForce.calculateForceAndEnergy(numStretchBends, posData, particle1, particle2, particle3,
lengthABParameters, lengthCBParameters, angleParameters, k1Parameters,
k2Parameters, forceData);
return static_cast<double>(energy);
......@@ -357,11 +357,11 @@ void ReferenceCalcAmoebaStretchBendForceKernel::copyParametersToContext(ContextI
force.getStretchBendParameters(i, particle1Index, particle2Index, particle3Index, lengthAB, lengthCB, angle, k1, k2);
if (particle1Index != particle1[i] || particle2Index != particle2[i] || particle3Index != particle3[i])
throw OpenMMException("updateParametersInContext: The set of particles in a stretch-bend has changed");
lengthABParameters[i] = (RealOpenMM) lengthAB;
lengthCBParameters[i] = (RealOpenMM) lengthCB;
angleParameters[i] = (RealOpenMM) angle;
k1Parameters[i] = (RealOpenMM) k1;
k2Parameters[i] = (RealOpenMM) k2;
lengthABParameters[i] = lengthAB;
lengthCBParameters[i] = lengthCB;
angleParameters[i] = angle;
k1Parameters[i] = k1;
k2Parameters[i] = k2;
}
}
......@@ -385,28 +385,28 @@ void ReferenceCalcAmoebaOutOfPlaneBendForceKernel::initialize(const System& syst
particle2.push_back(particle2Index);
particle3.push_back(particle3Index);
particle4.push_back(particle4Index);
kParameters.push_back(static_cast<RealOpenMM>(k));
kParameters.push_back(k);
}
globalOutOfPlaneBendAngleCubic = static_cast<RealOpenMM>(force.getAmoebaGlobalOutOfPlaneBendCubic());
globalOutOfPlaneBendAngleQuartic = static_cast<RealOpenMM>(force.getAmoebaGlobalOutOfPlaneBendQuartic());
globalOutOfPlaneBendAnglePentic = static_cast<RealOpenMM>(force.getAmoebaGlobalOutOfPlaneBendPentic());
globalOutOfPlaneBendAngleSextic = static_cast<RealOpenMM>(force.getAmoebaGlobalOutOfPlaneBendSextic());
globalOutOfPlaneBendAngleCubic = force.getAmoebaGlobalOutOfPlaneBendCubic();
globalOutOfPlaneBendAngleQuartic = force.getAmoebaGlobalOutOfPlaneBendQuartic();
globalOutOfPlaneBendAnglePentic = force.getAmoebaGlobalOutOfPlaneBendPentic();
globalOutOfPlaneBendAngleSextic = force.getAmoebaGlobalOutOfPlaneBendSextic();
usePeriodic = force.usesPeriodicBoundaryConditions();
}
double ReferenceCalcAmoebaOutOfPlaneBendForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
AmoebaReferenceOutOfPlaneBendForce amoebaReferenceOutOfPlaneBendForce;
if (usePeriodic)
amoebaReferenceOutOfPlaneBendForce.setPeriodic(extractBoxVectors(context));
RealOpenMM energy = amoebaReferenceOutOfPlaneBendForce.calculateForceAndEnergy(numOutOfPlaneBends, posData,
particle1, particle2, particle3, particle4,
kParameters,
globalOutOfPlaneBendAngleCubic,
globalOutOfPlaneBendAngleQuartic,
globalOutOfPlaneBendAnglePentic,
globalOutOfPlaneBendAngleSextic, forceData);
double energy = amoebaReferenceOutOfPlaneBendForce.calculateForceAndEnergy(numOutOfPlaneBends, posData,
particle1, particle2, particle3, particle4,
kParameters,
globalOutOfPlaneBendAngleCubic,
globalOutOfPlaneBendAngleQuartic,
globalOutOfPlaneBendAnglePentic,
globalOutOfPlaneBendAngleSextic, forceData);
return static_cast<double>(energy);
}
......@@ -422,7 +422,7 @@ void ReferenceCalcAmoebaOutOfPlaneBendForceKernel::copyParametersToContext(Conte
force.getOutOfPlaneBendParameters(i, particle1Index, particle2Index, particle3Index, particle4Index, k);
if (particle1Index != particle1[i] || particle2Index != particle2[i] || particle3Index != particle3[i] || particle4Index != particle4[i])
throw OpenMMException("updateParametersInContext: The set of particles in an out-of-plane bend has changed");
kParameters[i] = (RealOpenMM) k;
kParameters[i] = k;
}
}
......@@ -479,11 +479,11 @@ void ReferenceCalcAmoebaTorsionTorsionForceKernel::initialize(const System& syst
torsionTorsionGrids[ii][kk][jj].resize(grid[kk][jj].size());
if (reorder) {
for (unsigned int ll = 0; ll < grid[ll][jj].size(); ll++) {
torsionTorsionGrids[ii][kk][jj][ll] = static_cast<RealOpenMM>(reorderedGrid[kk][jj][ll]);
torsionTorsionGrids[ii][kk][jj][ll] = reorderedGrid[kk][jj][ll];
}
} else {
for (unsigned int ll = 0; ll < grid[ll][jj].size(); ll++) {
torsionTorsionGrids[ii][kk][jj][ll] = static_cast<RealOpenMM>(grid[kk][jj][ll]);
torsionTorsionGrids[ii][kk][jj][ll] = grid[kk][jj][ll];
}
}
}
......@@ -493,12 +493,12 @@ void ReferenceCalcAmoebaTorsionTorsionForceKernel::initialize(const System& syst
double ReferenceCalcAmoebaTorsionTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
AmoebaReferenceTorsionTorsionForce amoebaReferenceTorsionTorsionForce;
if (usePeriodic)
amoebaReferenceTorsionTorsionForce.setPeriodic(extractBoxVectors(context));
RealOpenMM energy = amoebaReferenceTorsionTorsionForce.calculateForceAndEnergy(numTorsionTorsions, posData,
double energy = amoebaReferenceTorsionTorsionForce.calculateForceAndEnergy(numTorsionTorsions, posData,
particle1, particle2, particle3, particle4, particle5,
chiralCheckAtom, gridIndices, torsionTorsionGrids, forceData);
return static_cast<double>(energy);
......@@ -535,7 +535,6 @@ void ReferenceCalcAmoebaMultipoleForceKernel::initialize(const System& system, c
int dipoleIndex = 0;
int quadrupoleIndex = 0;
int maxCovalentRange = 0;
double totalCharge = 0.0;
for (int ii = 0; ii < numMultipoles; ii++) {
......@@ -554,24 +553,24 @@ void ReferenceCalcAmoebaMultipoleForceKernel::initialize(const System& system, c
multipoleAtomXs[ii] = multipoleAtomX;
multipoleAtomYs[ii] = multipoleAtomY;
charges[ii] = static_cast<RealOpenMM>(charge);
tholes[ii] = static_cast<RealOpenMM>(tholeD);
dampingFactors[ii] = static_cast<RealOpenMM>(dampingFactorD);
polarity[ii] = static_cast<RealOpenMM>(polarityD);
charges[ii] = charge;
tholes[ii] = tholeD;
dampingFactors[ii] = dampingFactorD;
polarity[ii] = polarityD;
dipoles[dipoleIndex++] = static_cast<RealOpenMM>(dipolesD[0]);
dipoles[dipoleIndex++] = static_cast<RealOpenMM>(dipolesD[1]);
dipoles[dipoleIndex++] = static_cast<RealOpenMM>(dipolesD[2]);
dipoles[dipoleIndex++] = dipolesD[0];
dipoles[dipoleIndex++] = dipolesD[1];
dipoles[dipoleIndex++] = dipolesD[2];
quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[0]);
quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[1]);
quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[2]);
quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[3]);
quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[4]);
quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[5]);
quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[6]);
quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[7]);
quadrupoles[quadrupoleIndex++] = static_cast<RealOpenMM>(quadrupolesD[8]);
quadrupoles[quadrupoleIndex++] = quadrupolesD[0];
quadrupoles[quadrupoleIndex++] = quadrupolesD[1];
quadrupoles[quadrupoleIndex++] = quadrupolesD[2];
quadrupoles[quadrupoleIndex++] = quadrupolesD[3];
quadrupoles[quadrupoleIndex++] = quadrupolesD[4];
quadrupoles[quadrupoleIndex++] = quadrupolesD[5];
quadrupoles[quadrupoleIndex++] = quadrupolesD[6];
quadrupoles[quadrupoleIndex++] = quadrupolesD[7];
quadrupoles[quadrupoleIndex++] = quadrupolesD[8];
// covalent info
......@@ -646,7 +645,7 @@ AmoebaReferenceMultipoleForce* ReferenceCalcAmoebaMultipoleForceKernel::setupAmo
amoebaReferenceGeneralizedKirkwoodForce->setIncludeCavityTerm(gkKernel->getIncludeCavityTerm());
amoebaReferenceGeneralizedKirkwoodForce->setDirectPolarization(gkKernel->getDirectPolarization());
vector<RealOpenMM> parameters;
vector<double> parameters;
gkKernel->getAtomicRadii(parameters);
amoebaReferenceGeneralizedKirkwoodForce->setAtomicRadii(parameters);
......@@ -658,7 +657,7 @@ AmoebaReferenceMultipoleForce* ReferenceCalcAmoebaMultipoleForceKernel::setupAmo
// calculate Grycuk Born radii
vector<RealVec>& posData = extractPositions(context);
vector<Vec3>& posData = extractPositions(context);
amoebaReferenceGeneralizedKirkwoodForce->calculateGrycukBornRadii(posData);
amoebaReferenceMultipoleForce = new AmoebaReferenceGeneralizedKirkwoodMultipoleForce(amoebaReferenceGeneralizedKirkwoodForce);
......@@ -669,7 +668,7 @@ AmoebaReferenceMultipoleForce* ReferenceCalcAmoebaMultipoleForceKernel::setupAmo
amoebaReferencePmeMultipoleForce->setAlphaEwald(alphaEwald);
amoebaReferencePmeMultipoleForce->setCutoffDistance(cutoffDistance);
amoebaReferencePmeMultipoleForce->setPmeGridDimensions(pmeGridDimension);
RealVec* boxVectors = extractBoxVectors(context);
Vec3* boxVectors = extractBoxVectors(context);
double minAllowedSize = 1.999999*cutoffDistance;
if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize) {
throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
......@@ -704,12 +703,12 @@ double ReferenceCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bo
AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
RealOpenMM energy = amoebaReferenceMultipoleForce->calculateForceAndEnergy(posData, charges, dipoles, quadrupoles, tholes,
dampingFactors, polarity, axisTypes,
multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
multipoleAtomCovalentInfo, forceData);
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
double energy = amoebaReferenceMultipoleForce->calculateForceAndEnergy(posData, charges, dipoles, quadrupoles, tholes,
dampingFactors, polarity, axisTypes,
multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
multipoleAtomCovalentInfo, forceData);
delete amoebaReferenceMultipoleForce;
......@@ -723,11 +722,11 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getInducedDipoles(ContextImpl& con
// Create an AmoebaReferenceMultipoleForce to do the calculation.
AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
vector<RealVec>& posData = extractPositions(context);
vector<Vec3>& posData = extractPositions(context);
// Retrieve the induced dipoles.
vector<RealVec> inducedDipoles;
vector<Vec3> inducedDipoles;
amoebaReferenceMultipoleForce->calculateInducedDipoles(posData, charges, dipoles, quadrupoles, tholes,
dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, multipoleAtomCovalentInfo, inducedDipoles);
for (int i = 0; i < numParticles; i++)
......@@ -742,11 +741,11 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getLabFramePermanentDipoles(Contex
// Create an AmoebaReferenceMultipoleForce to do the calculation.
AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
vector<RealVec>& posData = extractPositions(context);
vector<Vec3>& posData = extractPositions(context);
// Retrieve the permanent dipoles in the lab frame.
vector<RealVec> labFramePermanentDipoles;
vector<Vec3> labFramePermanentDipoles;
amoebaReferenceMultipoleForce->calculateLabFramePermanentDipoles(posData, charges, dipoles, quadrupoles, tholes,
dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, multipoleAtomCovalentInfo, labFramePermanentDipoles);
for (int i = 0; i < numParticles; i++)
......@@ -762,11 +761,11 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getTotalDipoles(ContextImpl& conte
// Create an AmoebaReferenceMultipoleForce to do the calculation.
AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
vector<RealVec>& posData = extractPositions(context);
vector<Vec3>& posData = extractPositions(context);
// Retrieve the permanent dipoles in the lab frame.
vector<RealVec> totalDipoles;
vector<Vec3> totalDipoles;
amoebaReferenceMultipoleForce->calculateTotalDipoles(posData, charges, dipoles, quadrupoles, tholes,
dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, multipoleAtomCovalentInfo, totalDipoles);
......@@ -781,9 +780,9 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextI
std::vector< double >& outputElectrostaticPotential) {
AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
vector<RealVec>& posData = extractPositions(context);
vector<RealVec> grid(inputGrid.size());
vector<RealOpenMM> potential(inputGrid.size());
vector<Vec3>& posData = extractPositions(context);
vector<Vec3> grid(inputGrid.size());
vector<double> potential(inputGrid.size());
for (unsigned int ii = 0; ii < inputGrid.size(); ii++) {
grid[ii] = inputGrid[ii];
}
......@@ -805,13 +804,13 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextI
// retrieve masses
const System& system = context.getSystem();
vector<RealOpenMM> masses;
vector<double> masses;
for (int i = 0; i < system.getNumParticles(); ++i) {
masses.push_back(static_cast<RealOpenMM>(system.getParticleMass(i)));
masses.push_back(system.getParticleMass(i));
}
AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
vector<RealVec>& posData = extractPositions(context);
vector<Vec3>& posData = extractPositions(context);
amoebaReferenceMultipoleForce->calculateAmoebaSystemMultipoleMoments(masses, posData, charges, dipoles, quadrupoles, tholes,
dampingFactors, polarity, axisTypes,
multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
......@@ -838,22 +837,22 @@ void ReferenceCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImp
multipoleAtomZs[i] = multipoleAtomZ;
multipoleAtomXs[i] = multipoleAtomX;
multipoleAtomYs[i] = multipoleAtomY;
charges[i] = (RealOpenMM) charge;
tholes[i] = (RealOpenMM) tholeD;
dampingFactors[i] = (RealOpenMM) dampingFactorD;
polarity[i] = (RealOpenMM) polarityD;
dipoles[dipoleIndex++] = (RealOpenMM) dipolesD[0];
dipoles[dipoleIndex++] = (RealOpenMM) dipolesD[1];
dipoles[dipoleIndex++] = (RealOpenMM) dipolesD[2];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[0];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[1];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[2];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[3];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[4];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[5];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[6];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[7];
quadrupoles[quadrupoleIndex++] = (RealOpenMM) quadrupolesD[8];
charges[i] = charge;
tholes[i] = tholeD;
dampingFactors[i] = dampingFactorD;
polarity[i] = polarityD;
dipoles[dipoleIndex++] = dipolesD[0];
dipoles[dipoleIndex++] = dipolesD[1];
dipoles[dipoleIndex++] = dipolesD[2];
quadrupoles[quadrupoleIndex++] = quadrupolesD[0];
quadrupoles[quadrupoleIndex++] = quadrupolesD[1];
quadrupoles[quadrupoleIndex++] = quadrupolesD[2];
quadrupoles[quadrupoleIndex++] = quadrupolesD[3];
quadrupoles[quadrupoleIndex++] = quadrupolesD[4];
quadrupoles[quadrupoleIndex++] = quadrupolesD[5];
quadrupoles[quadrupoleIndex++] = quadrupolesD[6];
quadrupoles[quadrupoleIndex++] = quadrupolesD[7];
quadrupoles[quadrupoleIndex++] = quadrupolesD[8];
}
}
......@@ -889,37 +888,37 @@ int ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDirectPolarization() c
return directPolarization;
}
RealOpenMM ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSoluteDielectric() const {
double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSoluteDielectric() const {
return soluteDielectric;
}
RealOpenMM ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSolventDielectric() const {
double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSolventDielectric() const {
return solventDielectric;
}
RealOpenMM ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDielectricOffset() const {
double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDielectricOffset() const {
return dielectricOffset;
}
RealOpenMM ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getProbeRadius() const {
double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getProbeRadius() const {
return probeRadius;
}
RealOpenMM ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSurfaceAreaFactor() const {
double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSurfaceAreaFactor() const {
return surfaceAreaFactor;
}
void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getAtomicRadii(vector<RealOpenMM>& outputAtomicRadii) const {
void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getAtomicRadii(vector<double>& outputAtomicRadii) const {
outputAtomicRadii.resize(atomicRadii.size());
copy(atomicRadii.begin(), atomicRadii.end(), outputAtomicRadii.begin());
}
void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getScaleFactors(vector<RealOpenMM>& outputScaleFactors) const {
void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getScaleFactors(vector<double>& outputScaleFactors) const {
outputScaleFactors.resize(scaleFactors.size());
copy(scaleFactors.begin(), scaleFactors.end(), outputScaleFactors.begin());
}
void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getCharges(vector<RealOpenMM>& outputCharges) const {
void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getCharges(vector<double>& outputCharges) const {
outputCharges.resize(charges.size());
copy(charges.begin(), charges.end(), outputCharges.begin());
}
......@@ -947,9 +946,9 @@ void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System&
double particleCharge, particleRadius, scalingFactor;
force.getParticleParameters(ii, particleCharge, particleRadius, scalingFactor);
atomicRadii.push_back(static_cast<RealOpenMM>(particleRadius));
scaleFactors.push_back(static_cast<RealOpenMM>(scalingFactor));
charges.push_back(static_cast<RealOpenMM>(particleCharge));
atomicRadii.push_back(particleRadius);
scaleFactors.push_back(scalingFactor);
charges.push_back(particleCharge);
// Make sure the charge matches the one specified by the AmoebaMultipoleForce.
......@@ -963,11 +962,11 @@ void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System&
}
includeCavityTerm = force.getIncludeCavityTerm();
soluteDielectric = static_cast<RealOpenMM>(force.getSoluteDielectric());
solventDielectric = static_cast<RealOpenMM>(force.getSolventDielectric());
dielectricOffset = static_cast<RealOpenMM>(0.009);
probeRadius = static_cast<RealOpenMM>(force.getProbeRadius()),
surfaceAreaFactor = static_cast<RealOpenMM>(force.getSurfaceAreaFactor());
soluteDielectric = force.getSoluteDielectric();
solventDielectric = force.getSolventDielectric();
dielectricOffset = 0.009;
probeRadius = force.getProbeRadius(),
surfaceAreaFactor = force.getSurfaceAreaFactor();
directPolarization = amoebaMultipoleForce->getPolarizationType() == AmoebaMultipoleForce::Direct ? 1 : 0;
}
......@@ -1030,9 +1029,9 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
}
indexIVs[ii] = indexIV;
sigmas[ii] = static_cast<RealOpenMM>(sigma);
epsilons[ii] = static_cast<RealOpenMM>(epsilon);
reductions[ii] = static_cast<RealOpenMM>(reduction);
sigmas[ii] = sigma;
epsilons[ii] = epsilon;
reductions[ii] = reduction;
}
sigmaCombiningRule = force.getSigmaCombiningRule();
epsilonCombiningRule = force.getEpsilonCombiningRule();
......@@ -1046,16 +1045,16 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
AmoebaReferenceVdwForce vdwForce(sigmaCombiningRule, epsilonCombiningRule);
RealOpenMM energy;
double energy;
if (useCutoff) {
vdwForce.setCutoff(cutoff);
computeNeighborListVoxelHash(*neighborList, numParticles, posData, allExclusions, extractBoxVectors(context), usePBC, cutoff, 0.0);
if (usePBC) {
vdwForce.setNonbondedMethod(AmoebaReferenceVdwForce::CutoffPeriodic);
RealVec* boxVectors = extractBoxVectors(context);
Vec3* boxVectors = extractBoxVectors(context);
double minAllowedSize = 1.999999*cutoff;
if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize) {
throw OpenMMException("The periodic box size has decreased to less than twice the cutoff.");
......@@ -1084,9 +1083,9 @@ void ReferenceCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& con
double sigma, epsilon, reduction;
force.getParticleParameters(i, indexIV, sigma, epsilon, reduction);
indexIVs[i] = indexIV;
sigmas[i] = (RealOpenMM) sigma;
epsilons[i] = (RealOpenMM) epsilon;
reductions[i]= (RealOpenMM) reduction;
sigmas[i] = sigma;
epsilons[i] = epsilon;
reductions[i]= reduction;
}
}
......@@ -1113,27 +1112,27 @@ void ReferenceCalcAmoebaWcaDispersionForceKernel::initialize(const System& syste
double radius, epsilon;
force.getParticleParameters(ii, radius, epsilon);
radii[ii] = static_cast<RealOpenMM>(radius);
epsilons[ii] = static_cast<RealOpenMM>(epsilon);
radii[ii] = radius;
epsilons[ii] = epsilon;
}
totalMaximumDispersionEnergy = static_cast<RealOpenMM>(AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(force));
totalMaximumDispersionEnergy = AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(force);
epso = static_cast<RealOpenMM>(force.getEpso());
epsh = static_cast<RealOpenMM>(force.getEpsh());
rmino = static_cast<RealOpenMM>(force.getRmino());
rminh = static_cast<RealOpenMM>(force.getRminh());
awater = static_cast<RealOpenMM>(force.getAwater());
shctd = static_cast<RealOpenMM>(force.getShctd());
dispoff = static_cast<RealOpenMM>(force.getDispoff());
slevy = static_cast<RealOpenMM>(force.getSlevy());
epso = force.getEpso();
epsh = force.getEpsh();
rmino = force.getRmino();
rminh = force.getRminh();
awater = force.getAwater();
shctd = force.getShctd();
dispoff = force.getDispoff();
slevy = force.getSlevy();
}
double ReferenceCalcAmoebaWcaDispersionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& forceData = extractForces(context);
AmoebaReferenceWcaDispersionForce amoebaReferenceWcaDispersionForce(epso, epsh, rmino, rminh, awater, shctd, dispoff, slevy);
RealOpenMM energy = amoebaReferenceWcaDispersionForce.calculateForceAndEnergy(numParticles, posData, radii, epsilons, totalMaximumDispersionEnergy, forceData);
double energy = amoebaReferenceWcaDispersionForce.calculateForceAndEnergy(numParticles, posData, radii, epsilons, totalMaximumDispersionEnergy, forceData);
return static_cast<double>(energy);
}
......@@ -1146,8 +1145,8 @@ void ReferenceCalcAmoebaWcaDispersionForceKernel::copyParametersToContext(Contex
for (int i = 0; i < numParticles; ++i) {
double radius, epsilon;
force.getParticleParameters(i, radius, epsilon);
radii[i] = (RealOpenMM) radius;
epsilons[i] = (RealOpenMM) epsilon;
radii[i] = radius;
epsilons[i] = epsilon;
}
totalMaximumDispersionEnergy = (RealOpenMM) AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(force);
totalMaximumDispersionEnergy = AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(force);
}
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ a783b996
......@@ -72,10 +72,10 @@ private:
int numBonds;
std::vector<int> particle1;
std::vector<int> particle2;
std::vector<RealOpenMM> length;
std::vector<RealOpenMM> kQuadratic;
RealOpenMM globalBondCubic;
RealOpenMM globalBondQuartic;
std::vector<double> length;
std::vector<double> kQuadratic;
double globalBondCubic;
double globalBondQuartic;
const System& system;
bool usePeriodic;
};
......@@ -115,12 +115,12 @@ private:
std::vector<int> particle1;
std::vector<int> particle2;
std::vector<int> particle3;
std::vector<RealOpenMM> angle;
std::vector<RealOpenMM> kQuadratic;
RealOpenMM globalAngleCubic;
RealOpenMM globalAngleQuartic;
RealOpenMM globalAnglePentic;
RealOpenMM globalAngleSextic;
std::vector<double> angle;
std::vector<double> kQuadratic;
double globalAngleCubic;
double globalAngleQuartic;
double globalAnglePentic;
double globalAngleSextic;
const System& system;
bool usePeriodic;
};
......@@ -161,12 +161,12 @@ private:
std::vector<int> particle2;
std::vector<int> particle3;
std::vector<int> particle4;
std::vector<RealOpenMM> angle;
std::vector<RealOpenMM> kQuadratic;
RealOpenMM globalInPlaneAngleCubic;
RealOpenMM globalInPlaneAngleQuartic;
RealOpenMM globalInPlaneAnglePentic;
RealOpenMM globalInPlaneAngleSextic;
std::vector<double> angle;
std::vector<double> kQuadratic;
double globalInPlaneAngleCubic;
double globalInPlaneAngleQuartic;
double globalInPlaneAnglePentic;
double globalInPlaneAngleSextic;
const System& system;
bool usePeriodic;
};
......@@ -209,7 +209,7 @@ private:
std::vector<int> particle4;
std::vector<int> particle5;
std::vector<int> particle6;
std::vector<RealOpenMM> kTorsion;
std::vector<double> kTorsion;
const System& system;
bool usePeriodic;
};
......@@ -249,11 +249,11 @@ private:
std::vector<int> particle1;
std::vector<int> particle2;
std::vector<int> particle3;
std::vector<RealOpenMM> lengthABParameters;
std::vector<RealOpenMM> lengthCBParameters;
std::vector<RealOpenMM> angleParameters;
std::vector<RealOpenMM> k1Parameters;
std::vector<RealOpenMM> k2Parameters;
std::vector<double> lengthABParameters;
std::vector<double> lengthCBParameters;
std::vector<double> angleParameters;
std::vector<double> k1Parameters;
std::vector<double> k2Parameters;
const System& system;
bool usePeriodic;
};
......@@ -294,11 +294,11 @@ private:
std::vector<int> particle2;
std::vector<int> particle3;
std::vector<int> particle4;
std::vector<RealOpenMM> kParameters;
RealOpenMM globalOutOfPlaneBendAngleCubic;
RealOpenMM globalOutOfPlaneBendAngleQuartic;
RealOpenMM globalOutOfPlaneBendAnglePentic;
RealOpenMM globalOutOfPlaneBendAngleSextic;
std::vector<double> kParameters;
double globalOutOfPlaneBendAngleCubic;
double globalOutOfPlaneBendAngleQuartic;
double globalOutOfPlaneBendAnglePentic;
double globalOutOfPlaneBendAngleSextic;
const System& system;
bool usePeriodic;
};
......@@ -337,7 +337,7 @@ private:
std::vector<int> gridIndices;
int numTorsionTorsionGrids;
std::vector< std::vector< std::vector< std::vector<RealOpenMM> > > > torsionTorsionGrids;
std::vector< std::vector< std::vector< std::vector<double> > > > torsionTorsionGrids;
const System& system;
bool usePeriodic;
......@@ -439,12 +439,12 @@ private:
int numMultipoles;
AmoebaMultipoleForce::NonbondedMethod nonbondedMethod;
AmoebaMultipoleForce::PolarizationType polarizationType;
std::vector<RealOpenMM> charges;
std::vector<RealOpenMM> dipoles;
std::vector<RealOpenMM> quadrupoles;
std::vector<RealOpenMM> tholes;
std::vector<RealOpenMM> dampingFactors;
std::vector<RealOpenMM> polarity;
std::vector<double> charges;
std::vector<double> dipoles;
std::vector<double> quadrupoles;
std::vector<double> tholes;
std::vector<double> dampingFactors;
std::vector<double> polarity;
std::vector<int> axisTypes;
std::vector<int> multipoleAtomZs;
std::vector<int> multipoleAtomXs;
......@@ -452,12 +452,12 @@ private:
std::vector< std::vector< std::vector<int> > > multipoleAtomCovalentInfo;
int mutualInducedMaxIterations;
RealOpenMM mutualInducedTargetEpsilon;
double mutualInducedTargetEpsilon;
std::vector<double> extrapolationCoefficients;
bool usePme;
RealOpenMM alphaEwald;
RealOpenMM cutoffDistance;
double alphaEwald;
double cutoffDistance;
std::vector<int> pmeGridDimension;
const System& system;
......@@ -501,9 +501,9 @@ private:
double dispersionCoefficient;
std::vector<int> indexIVs;
std::vector< std::set<int> > allExclusions;
std::vector<RealOpenMM> sigmas;
std::vector<RealOpenMM> epsilons;
std::vector<RealOpenMM> reductions;
std::vector<double> sigmas;
std::vector<double> epsilons;
std::vector<double> reductions;
std::string sigmaCombiningRule;
std::string epsilonCombiningRule;
const System& system;
......@@ -543,17 +543,17 @@ public:
private:
int numParticles;
std::vector<RealOpenMM> radii;
std::vector<RealOpenMM> epsilons;
RealOpenMM epso;
RealOpenMM epsh;
RealOpenMM rmino;
RealOpenMM rminh;
RealOpenMM awater;
RealOpenMM shctd;
RealOpenMM dispoff;
RealOpenMM slevy;
RealOpenMM totalMaximumDispersionEnergy;
std::vector<double> radii;
std::vector<double> epsilons;
double epso;
double epsh;
double rmino;
double rminh;
double awater;
double shctd;
double dispoff;
double slevy;
double totalMaximumDispersionEnergy;
const System& system;
};
......@@ -609,7 +609,7 @@ public:
* @return soluteDielectric
*
*/
RealOpenMM getSoluteDielectric() const;
double getSoluteDielectric() const;
/**
* Get the solvent dielectric.
......@@ -617,7 +617,7 @@ public:
* @return solventDielectric
*
*/
RealOpenMM getSolventDielectric() const;
double getSolventDielectric() const;
/**
* Get the dielectric offset.
......@@ -625,7 +625,7 @@ public:
* @return dielectricOffset
*
*/
RealOpenMM getDielectricOffset() const;
double getDielectricOffset() const;
/**
* Get the probe radius.
......@@ -633,7 +633,7 @@ public:
* @return probeRadius
*
*/
RealOpenMM getProbeRadius() const;
double getProbeRadius() const;
/**
* Get the surface area factor.
......@@ -641,7 +641,7 @@ public:
* @return surfaceAreaFactor
*
*/
RealOpenMM getSurfaceAreaFactor() const;
double getSurfaceAreaFactor() const;
/**
* Get the vector of particle radii.
......@@ -649,7 +649,7 @@ public:
* @param atomicRadii vector of atomic radii
*
*/
void getAtomicRadii(std::vector<RealOpenMM>& atomicRadii) const;
void getAtomicRadii(std::vector<double>& atomicRadii) const;
/**
* Get the vector of scale factors.
......@@ -657,7 +657,7 @@ public:
* @param scaleFactors vector of scale factors
*
*/
void getScaleFactors(std::vector<RealOpenMM>& scaleFactors) const;
void getScaleFactors(std::vector<double>& scaleFactors) const;
/**
* Get the vector of charges.
......@@ -665,7 +665,7 @@ public:
* @param charges vector of charges
*
*/
void getCharges(std::vector<RealOpenMM>& charges) const;
void getCharges(std::vector<double>& charges) const;
/**
* Copy changed parameters over to a context.
......@@ -678,14 +678,14 @@ public:
private:
int numParticles;
std::vector<RealOpenMM> atomicRadii;
std::vector<RealOpenMM> scaleFactors;
std::vector<RealOpenMM> charges;
RealOpenMM soluteDielectric;
RealOpenMM solventDielectric;
RealOpenMM dielectricOffset;
RealOpenMM probeRadius;
RealOpenMM surfaceAreaFactor;
std::vector<double> atomicRadii;
std::vector<double> scaleFactors;
std::vector<double> charges;
double soluteDielectric;
double solventDielectric;
double dielectricOffset;
double probeRadius;
double surfaceAreaFactor;
int includeCavityTerm;
int directPolarization;
const System& system;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceAngleForce.cpp
View file @ a783b996
......@@ -24,11 +24,12 @@
#include "AmoebaReferenceForce.h"
#include "AmoebaReferenceAngleForce.h"
#include "SimTKOpenMMRealType.h"
using std::vector;
using namespace OpenMM;
void AmoebaReferenceAngleForce::setPeriodic(OpenMM::RealVec* vectors) {
void AmoebaReferenceAngleForce::setPeriodic(OpenMM::Vec3* vectors) {
usePeriodic = true;
boxVectors[0] = vectors[0];
boxVectors[1] = vectors[1];
......@@ -53,49 +54,35 @@ void AmoebaReferenceAngleForce::setPeriodic(OpenMM::RealVec* vectors) {
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceAngleForce::getPrefactorsGivenAngleCosine(RealOpenMM cosine,
RealOpenMM idealAngle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM* dEdR) const {
double AmoebaReferenceAngleForce::getPrefactorsGivenAngleCosine(double cosine,
double idealAngle, double angleK,
double angleCubic, double angleQuartic,
double anglePentic, double angleSextic,
double* dEdR) const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
static const RealOpenMM four = 4.0;
static const RealOpenMM five = 5.0;
static const RealOpenMM six = 6.0;
// static const std::string methodName = "AmoebaReferenceAngleForce::getPrefactorsGivenAngleCosine";
// ---------------------------------------------------------------------------------------
RealOpenMM angle;
if (cosine >= one) {
angle = zero;
} else if (cosine <= -one) {
double angle;
if (cosine >= 1.0) {
angle = 0.0;
} else if (cosine <= -1.0) {
angle = RADIAN*PI_M;
} else {
angle = RADIAN*ACOS(cosine);
}
RealOpenMM deltaIdeal = angle - idealAngle;
RealOpenMM deltaIdeal2 = deltaIdeal*deltaIdeal;
RealOpenMM deltaIdeal3 = deltaIdeal*deltaIdeal2;
RealOpenMM deltaIdeal4 = deltaIdeal2*deltaIdeal2;
double deltaIdeal = angle - idealAngle;
double deltaIdeal2 = deltaIdeal*deltaIdeal;
double deltaIdeal3 = deltaIdeal*deltaIdeal2;
double deltaIdeal4 = deltaIdeal2*deltaIdeal2;
*dEdR = (two + three*angleCubic*deltaIdeal +
four*angleQuartic*deltaIdeal2 +
five*anglePentic*deltaIdeal3 +
six*angleSextic*deltaIdeal4);
*dEdR = (2.0 + 3.0*angleCubic*deltaIdeal +
4.0*angleQuartic*deltaIdeal2 +
5.0*anglePentic*deltaIdeal3 +
6.0*angleSextic*deltaIdeal4);
*dEdR *= RADIAN*angleK*deltaIdeal;
*dEdR *= RADIAN*angleK*deltaIdeal;
RealOpenMM energy = 1.0f + angleCubic*deltaIdeal + angleQuartic*deltaIdeal2 +
anglePentic*deltaIdeal3 + angleSextic*deltaIdeal4;
energy *= angleK*deltaIdeal2;
double energy = 1.0f + angleCubic*deltaIdeal + angleQuartic*deltaIdeal2 +
anglePentic*deltaIdeal3 + angleSextic*deltaIdeal4;
energy *= angleK*deltaIdeal2;
return energy;
......@@ -120,60 +107,49 @@ RealOpenMM AmoebaReferenceAngleForce::getPrefactorsGivenAngleCosine(RealOpenMM c
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceAngleForce::calculateAngleIxn(const RealVec& positionAtomA, const RealVec& positionAtomB,
const RealVec& positionAtomC,
RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealVec* forces) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceAngleForce::calculateAngleIxn";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM onePt5 = 1.5;
static const RealOpenMM two = 2.0;
// ---------------------------------------------------------------------------------------
double AmoebaReferenceAngleForce::calculateAngleIxn(const Vec3& positionAtomA, const Vec3& positionAtomB,
const Vec3& positionAtomC,
double angle, double angleK,
double angleCubic, double angleQuartic,
double anglePentic, double angleSextic,
Vec3* forces) const {
std::vector<RealOpenMM> deltaR[2];
std::vector<double> deltaR[2];
if (usePeriodic)
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomA, positionAtomB, deltaR[0], boxVectors);
else
AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomB, deltaR[0]);
RealOpenMM rAB2 = AmoebaReferenceForce::getNormSquared3(deltaR[0]);
RealOpenMM rAB = SQRT(rAB2);
double rAB2 = AmoebaReferenceForce::getNormSquared3(deltaR[0]);
double rAB = SQRT(rAB2);
if (usePeriodic)
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomC, positionAtomB, deltaR[1], boxVectors);
else
AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomB, deltaR[1]);
RealOpenMM rCB2 = AmoebaReferenceForce::getNormSquared3(deltaR[1]);
RealOpenMM rCB = SQRT(rCB2);
double rCB2 = AmoebaReferenceForce::getNormSquared3(deltaR[1]);
double rCB = SQRT(rCB2);
if (rAB <= zero || rCB <= zero) {
return zero;
if (rAB <= 0.0 || rCB <= 0.0) {
return 0.0;
}
std::vector<RealOpenMM> pVector(3);
std::vector<double> pVector(3);
AmoebaReferenceForce::getCrossProduct(deltaR[0], deltaR[1], pVector);
RealOpenMM rp = AmoebaReferenceForce::getNorm3(pVector);
double rp = AmoebaReferenceForce::getNorm3(pVector);
if (rp < 1.0e-06) {
rp = 1.0e-06;
}
RealOpenMM dot = AmoebaReferenceForce::getDotProduct3(deltaR[0], deltaR[1]);
RealOpenMM cosine = dot/(rAB*rCB);
double dot = AmoebaReferenceForce::getDotProduct3(deltaR[0], deltaR[1]);
double cosine = dot/(rAB*rCB);
RealOpenMM dEdR;
RealOpenMM energy = getPrefactorsGivenAngleCosine(cosine, angle, angleK, angleCubic, angleQuartic,
double dEdR;
double energy = getPrefactorsGivenAngleCosine(cosine, angle, angleK, angleCubic, angleQuartic,
anglePentic, angleSextic, &dEdR);
RealOpenMM termA = dEdR/(rAB2*rp);
RealOpenMM termC = -dEdR/(rCB2*rp);
double termA = dEdR/(rAB2*rp);
double termC = -dEdR/(rCB2*rp);
std::vector<RealOpenMM> deltaCrossP[3];
std::vector<double> deltaCrossP[3];
deltaCrossP[0].resize(3);
deltaCrossP[1].resize(3);
deltaCrossP[2].resize(3);
......@@ -196,27 +172,27 @@ RealOpenMM AmoebaReferenceAngleForce::calculateAngleIxn(const RealVec& positionA
return energy;
}
RealOpenMM AmoebaReferenceAngleForce::calculateForceAndEnergy(int numAngles, vector<RealVec>& posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<RealOpenMM>& angle,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM angleCubic,
RealOpenMM angleQuartic,
RealOpenMM anglePentic,
RealOpenMM angleSextic,
vector<RealVec>& forceData) const {
RealOpenMM energy = 0.0;
double AmoebaReferenceAngleForce::calculateForceAndEnergy(int numAngles, vector<Vec3>& posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<double>& angle,
const std::vector<double>& kQuadratic,
double angleCubic,
double angleQuartic,
double anglePentic,
double angleSextic,
vector<Vec3>& forceData) const {
double energy = 0.0;
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numAngles); ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
RealOpenMM idealAngle = angle[ii];
RealOpenMM angleK = kQuadratic[ii];
RealVec forces[3];
energy += calculateAngleIxn(posData[particle1Index], posData[particle2Index], posData[particle3Index],
idealAngle, angleK, angleCubic, angleQuartic, anglePentic, angleSextic, forces);
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
double idealAngle = angle[ii];
double angleK = kQuadratic[ii];
Vec3 forces[3];
energy += calculateAngleIxn(posData[particle1Index], posData[particle2Index], posData[particle3Index],
idealAngle, angleK, angleCubic, angleQuartic, anglePentic, angleSextic, forces);
for (unsigned int jj = 0; jj < 3; jj++) {
forceData[particle1Index][jj] += forces[0][jj];
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceAngleForce.h
View file @ a783b996
......@@ -25,7 +25,7 @@
#ifndef __AmoebaReferenceAngleForce_H__
#define __AmoebaReferenceAngleForce_H__
#include "RealVec.h"
#include "openmm/Vec3.h"
#include <vector>
namespace OpenMM {
......@@ -58,7 +58,7 @@ public:
--------------------------------------------------------------------------------------- */
void setPeriodic(OpenMM::RealVec* vectors);
void setPeriodic(OpenMM::Vec3* vectors);
/**---------------------------------------------------------------------------------------
......@@ -82,22 +82,22 @@ public:
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy(int numAngles, std::vector<OpenMM::RealVec>& posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<RealOpenMM>& angle,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM globalAngleCubic,
RealOpenMM globalAngleQuartic,
RealOpenMM globalAnglePentic,
RealOpenMM globalAngleSextic,
std::vector<OpenMM::RealVec>& forceData) const;
double calculateForceAndEnergy(int numAngles, std::vector<OpenMM::Vec3>& posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<double>& angle,
const std::vector<double>& kQuadratic,
double globalAngleCubic,
double globalAngleQuartic,
double globalAnglePentic,
double globalAngleSextic,
std::vector<OpenMM::Vec3>& forceData) const;
private:
bool usePeriodic;
RealVec boxVectors[3];
Vec3 boxVectors[3];
/**---------------------------------------------------------------------------------------
......@@ -117,10 +117,10 @@ private:
--------------------------------------------------------------------------------------- */
RealOpenMM getPrefactorsGivenAngleCosine(RealOpenMM cosine, RealOpenMM idealAngle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
RealOpenMM* dEdR) const;
double getPrefactorsGivenAngleCosine(double cosine, double idealAngle, double angleK,
double angleCubic, double angleQuartic,
double anglePentic, double angleSextic,
double* dEdR) const;
/**---------------------------------------------------------------------------------------
......@@ -141,12 +141,12 @@ private:
--------------------------------------------------------------------------------------- */
RealOpenMM calculateAngleIxn(const OpenMM::RealVec& positionAtomA, const OpenMM::RealVec& positionAtomB,
const OpenMM::RealVec& positionAtomC,
RealOpenMM angle, RealOpenMM angleK,
RealOpenMM angleCubic, RealOpenMM angleQuartic,
RealOpenMM anglePentic, RealOpenMM angleSextic,
OpenMM::RealVec* forces) const;
double calculateAngleIxn(const OpenMM::Vec3& positionAtomA, const OpenMM::Vec3& positionAtomB,
const OpenMM::Vec3& positionAtomC,
double angle, double angleK,
double angleCubic, double angleQuartic,
double anglePentic, double angleSextic,
OpenMM::Vec3* forces) const;
};
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceBondForce.cpp
View file @ a783b996
......@@ -28,7 +28,7 @@
using std::vector;
using namespace OpenMM;
void AmoebaReferenceBondForce::setPeriodic(OpenMM::RealVec* vectors) {
void AmoebaReferenceBondForce::setPeriodic(OpenMM::Vec3* vectors) {
usePeriodic = true;
boxVectors[0] = vectors[0];
boxVectors[1] = vectors[1];
......@@ -51,73 +51,61 @@ void AmoebaReferenceBondForce::setPeriodic(OpenMM::RealVec* vectors) {
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceBondForce::calculateBondIxn(const RealVec& positionAtomA, const RealVec& positionAtomB,
RealOpenMM bondLength, RealOpenMM bondK,
RealOpenMM bondCubic, RealOpenMM bondQuartic,
RealVec* forces) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceBondForce::calculateBondIxn";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM onePt5 = 1.5;
static const RealOpenMM two = 2.0;
// ---------------------------------------------------------------------------------------
double AmoebaReferenceBondForce::calculateBondIxn(const Vec3& positionAtomA, const Vec3& positionAtomB,
double bondLength, double bondK,
double bondCubic, double bondQuartic,
Vec3* forces) const {
// get deltaR, R2, and R between 2 atoms
std::vector<RealOpenMM> deltaR;
std::vector<double> deltaR;
if (usePeriodic)
AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomA, positionAtomB, deltaR, boxVectors);
else
AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomB, deltaR);
RealOpenMM r = AmoebaReferenceForce::getNorm3(deltaR);
double r = AmoebaReferenceForce::getNorm3(deltaR);
// deltaIdeal = r - r_0
RealOpenMM deltaIdeal = r - bondLength;
RealOpenMM deltaIdeal2 = deltaIdeal*deltaIdeal;
double deltaIdeal = r - bondLength;
double deltaIdeal2 = deltaIdeal*deltaIdeal;
RealOpenMM dEdR = (one + onePt5*bondCubic*deltaIdeal + two*bondQuartic*deltaIdeal2);
dEdR *= two*bondK*deltaIdeal;
dEdR = r > zero ? (dEdR/r) : zero;
double dEdR = (1.0 + 1.5*bondCubic*deltaIdeal + 2.0*bondQuartic*deltaIdeal2);
dEdR *= 2.0*bondK*deltaIdeal;
dEdR = r > 0.0 ? (dEdR/r) : 0.0;
forces[0][0] = dEdR*deltaR[0];
forces[0][1] = dEdR*deltaR[1];
forces[0][2] = dEdR*deltaR[2];
forces[0][0] = dEdR*deltaR[0];
forces[0][1] = dEdR*deltaR[1];
forces[0][2] = dEdR*deltaR[2];
dEdR *= -1.0;
forces[1][0] = dEdR*deltaR[0];
forces[1][1] = dEdR*deltaR[1];
forces[1][2] = dEdR*deltaR[2];
dEdR *= -1.0;
forces[1][0] = dEdR*deltaR[0];
forces[1][1] = dEdR*deltaR[1];
forces[1][2] = dEdR*deltaR[2];
RealOpenMM energy = bondK*deltaIdeal2*(one + bondCubic*deltaIdeal + bondQuartic*deltaIdeal2);
double energy = bondK*deltaIdeal2*(1.0 + bondCubic*deltaIdeal + bondQuartic*deltaIdeal2);
return energy;
}
RealOpenMM AmoebaReferenceBondForce::calculateForceAndEnergy(int numBonds,
vector<RealVec>& particlePositions,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<RealOpenMM>& length,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM globalBondCubic,
RealOpenMM globalBondQuartic,
vector<RealVec>& forceData) const {
RealOpenMM energy = 0.0;
double AmoebaReferenceBondForce::calculateForceAndEnergy(int numBonds,
vector<Vec3>& particlePositions,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<double>& length,
const std::vector<double>& kQuadratic,
double globalBondCubic,
double globalBondQuartic,
vector<Vec3>& forceData) const {
double energy = 0.0;
for (int ii = 0; ii < numBonds; ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
RealOpenMM bondLength = length[ii];
RealOpenMM bondK = kQuadratic[ii];
RealVec forces[2];
energy += calculateBondIxn(particlePositions[particle1Index], particlePositions[particle2Index],
bondLength, bondK, globalBondCubic, globalBondQuartic,
forces);
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
double bondLength = length[ii];
double bondK = kQuadratic[ii];
Vec3 forces[2];
energy += calculateBondIxn(particlePositions[particle1Index], particlePositions[particle2Index],
bondLength, bondK, globalBondCubic, globalBondQuartic,
forces);
for (int jj = 0; jj < 3; jj++) {
forceData[particle1Index][jj] += forces[0][jj];
......
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