Skip to content

GitLab

Commit 8469621f authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #1747 from peastman/realtype 

Eliminated RealOpenMM
parents b84e22ba 6813ca57
Hide whitespace changes
Inline Side-by-side
Showing with 433 additions and 491 deletions
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp
View file @ 8469621f
...@@ -27,6 +27,7 @@ ...@@ -27,6 +27,7 @@
#include "ReferenceForce.h" #include "ReferenceForce.h"
#include <algorithm> #include <algorithm>
#include <cctype> #include <cctype>
#include <cmath>
using std::vector; using std::vector;
using namespace OpenMM; using namespace OpenMM;
...@@ -76,7 +77,7 @@ double AmoebaReferenceVdwForce::getCutoff() const { ...@@ -76,7 +77,7 @@ double AmoebaReferenceVdwForce::getCutoff() const {
return _cutoff; return _cutoff;
} }
void AmoebaReferenceVdwForce::setPeriodicBox(OpenMM::RealVec* vectors) { void AmoebaReferenceVdwForce::setPeriodicBox(OpenMM::Vec3* vectors) {
_periodicBoxVectors[0] = vectors[0]; _periodicBoxVectors[0] = vectors[0];
_periodicBoxVectors[1] = vectors[1]; _periodicBoxVectors[1] = vectors[1];
_periodicBoxVectors[2] = vectors[2]; _periodicBoxVectors[2] = vectors[2];
...@@ -102,22 +103,19 @@ std::string AmoebaReferenceVdwForce::getSigmaCombiningRule() const { ...@@ -102,22 +103,19 @@ std::string AmoebaReferenceVdwForce::getSigmaCombiningRule() const {
return _sigmaCombiningRule; return _sigmaCombiningRule;
} }
RealOpenMM AmoebaReferenceVdwForce::arithmeticSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const { double AmoebaReferenceVdwForce::arithmeticSigmaCombiningRule(double sigmaI, double sigmaJ) const {
return (sigmaI + sigmaJ); return (sigmaI + sigmaJ);
} }
RealOpenMM AmoebaReferenceVdwForce::geometricSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const { double AmoebaReferenceVdwForce::geometricSigmaCombiningRule(double sigmaI, double sigmaJ) const {
return 2.0*SQRT(sigmaI*sigmaJ); return 2.0*sqrt(sigmaI*sigmaJ);
} }
RealOpenMM AmoebaReferenceVdwForce::cubicMeanSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const { double AmoebaReferenceVdwForce::cubicMeanSigmaCombiningRule(double sigmaI, double sigmaJ) const {
double sigmaI2 = sigmaI*sigmaI;
double sigmaJ2 = sigmaJ*sigmaJ;
const RealOpenMM zero = 0.0; return sigmaI != 0.0 && sigmaJ != 0.0 ? 2.0*(sigmaI2*sigmaI + sigmaJ2*sigmaJ)/(sigmaI2 + sigmaJ2) : 0.0;
RealOpenMM sigmaI2 = sigmaI*sigmaI;
RealOpenMM sigmaJ2 = sigmaJ*sigmaJ;
return sigmaI != zero && sigmaJ != 0.0 ? 2.0*(sigmaI2*sigmaI + sigmaJ2*sigmaJ)/(sigmaI2 + sigmaJ2) : zero;
} }
void AmoebaReferenceVdwForce::setEpsilonCombiningRule(const std::string& epsilonCombiningRule) { void AmoebaReferenceVdwForce::setEpsilonCombiningRule(const std::string& epsilonCombiningRule) {
...@@ -142,97 +140,83 @@ std::string AmoebaReferenceVdwForce::getEpsilonCombiningRule() const { ...@@ -142,97 +140,83 @@ std::string AmoebaReferenceVdwForce::getEpsilonCombiningRule() const {
return _epsilonCombiningRule; return _epsilonCombiningRule;
} }
RealOpenMM AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const { double AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule(double epsilonI, double epsilonJ) const {
return 0.5*(epsilonI + epsilonJ); return 0.5*(epsilonI + epsilonJ);
} }
RealOpenMM AmoebaReferenceVdwForce::geometricEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const { double AmoebaReferenceVdwForce::geometricEpsilonCombiningRule(double epsilonI, double epsilonJ) const {
return SQRT(epsilonI*epsilonJ); return sqrt(epsilonI*epsilonJ);
} }
RealOpenMM AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const { double AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ) const {
return (epsilonI != 0.0 && epsilonJ != 0.0) ? 2.0*(epsilonI*epsilonJ)/(epsilonI + epsilonJ) : 0.0; return (epsilonI != 0.0 && epsilonJ != 0.0) ? 2.0*(epsilonI*epsilonJ)/(epsilonI + epsilonJ) : 0.0;
} }
RealOpenMM AmoebaReferenceVdwForce::hhgEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const { double AmoebaReferenceVdwForce::hhgEpsilonCombiningRule(double epsilonI, double epsilonJ) const {
RealOpenMM denominator = SQRT(epsilonI) + SQRT(epsilonJ); double denominator = sqrt(epsilonI) + sqrt(epsilonJ);
return (epsilonI != 0.0 && epsilonJ != 0.0) ? 4.0*(epsilonI*epsilonJ)/(denominator*denominator) : 0.0; return (epsilonI != 0.0 && epsilonJ != 0.0) ? 4.0*(epsilonI*epsilonJ)/(denominator*denominator) : 0.0;
} }
void AmoebaReferenceVdwForce::addReducedForce(unsigned int particleI, unsigned int particleIV, void AmoebaReferenceVdwForce::addReducedForce(unsigned int particleI, unsigned int particleIV,
RealOpenMM reduction, RealOpenMM sign, double reduction, double sign,
Vec3& force, vector<RealVec>& forces) const { Vec3& force, vector<Vec3>& forces) const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
forces[particleI][0] += sign*force[0]*reduction; forces[particleI][0] += sign*force[0]*reduction;
forces[particleI][1] += sign*force[1]*reduction; forces[particleI][1] += sign*force[1]*reduction;
forces[particleI][2] += sign*force[2]*reduction; forces[particleI][2] += sign*force[2]*reduction;
forces[particleIV][0] += sign*force[0]*(one - reduction); forces[particleIV][0] += sign*force[0]*(1.0 - reduction);
forces[particleIV][1] += sign*force[1]*(one - reduction); forces[particleIV][1] += sign*force[1]*(1.0 - reduction);
forces[particleIV][2] += sign*force[2]*(one - reduction); forces[particleIV][2] += sign*force[2]*(1.0 - reduction);
} }
RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn(RealOpenMM combinedSigma, RealOpenMM combinedEpsilon, double AmoebaReferenceVdwForce::calculatePairIxn(double combinedSigma, double combinedEpsilon,
const Vec3& particleIPosition, const Vec3& particleIPosition,
const Vec3& particleJPosition, const Vec3& particleJPosition,
Vec3& force) const { Vec3& force) const {
// --------------------------------------------------------------------------------------- static const double dhal = 0.07;
static const double ghal = 0.12;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM seven = 7.0;
static const RealOpenMM dhal = 0.07;
static const RealOpenMM ghal = 0.12;
// ---------------------------------------------------------------------------------------
// get deltaR, R2, and R between 2 atoms // get deltaR, R2, and R between 2 atoms
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex]; double deltaR[ReferenceForce::LastDeltaRIndex];
if (_nonbondedMethod == CutoffPeriodic) if (_nonbondedMethod == CutoffPeriodic)
ReferenceForce::getDeltaRPeriodic(particleJPosition, particleIPosition, _periodicBoxVectors, deltaR); ReferenceForce::getDeltaRPeriodic(particleJPosition, particleIPosition, _periodicBoxVectors, deltaR);
else else
ReferenceForce::getDeltaR(particleJPosition, particleIPosition, deltaR); ReferenceForce::getDeltaR(particleJPosition, particleIPosition, deltaR);
RealOpenMM r_ij_2 = deltaR[ReferenceForce::R2Index]; double r_ij_2 = deltaR[ReferenceForce::R2Index];
RealOpenMM r_ij = deltaR[ReferenceForce::RIndex]; double r_ij = deltaR[ReferenceForce::RIndex];
RealOpenMM sigma_7 = combinedSigma*combinedSigma*combinedSigma; double sigma_7 = combinedSigma*combinedSigma*combinedSigma;
sigma_7 = sigma_7*sigma_7*combinedSigma; sigma_7 = sigma_7*sigma_7*combinedSigma;
RealOpenMM r_ij_6 = r_ij_2*r_ij_2*r_ij_2; double r_ij_6 = r_ij_2*r_ij_2*r_ij_2;
RealOpenMM r_ij_7 = r_ij_6*r_ij; double r_ij_7 = r_ij_6*r_ij;
RealOpenMM rho = r_ij_7 + ghal*sigma_7; double rho = r_ij_7 + ghal*sigma_7;
RealOpenMM tau = (dhal + one)/(r_ij + dhal*combinedSigma); double tau = (dhal + 1.0)/(r_ij + dhal*combinedSigma);
RealOpenMM tau_7 = tau*tau*tau; double tau_7 = tau*tau*tau;
tau_7 = tau_7*tau_7*tau; tau_7 = tau_7*tau_7*tau;
RealOpenMM dtau = tau/(dhal + one); double dtau = tau/(dhal + 1.0);
RealOpenMM ratio = (sigma_7/rho); double ratio = (sigma_7/rho);
RealOpenMM gtau = combinedEpsilon*tau_7*r_ij_6*(ghal+one)*ratio*ratio; double gtau = combinedEpsilon*tau_7*r_ij_6*(ghal+1.0)*ratio*ratio;
RealOpenMM energy = combinedEpsilon*tau_7*sigma_7*((ghal+one)*sigma_7/rho - two); double energy = combinedEpsilon*tau_7*sigma_7*((ghal+1.0)*sigma_7/rho - 2.0);
RealOpenMM dEdR = -seven*(dtau*energy + gtau); double dEdR = -7.0*(dtau*energy + gtau);
// tapering // tapering
if ((_nonbondedMethod == CutoffNonPeriodic || _nonbondedMethod == CutoffPeriodic) && r_ij > _taperCutoff) { if ((_nonbondedMethod == CutoffNonPeriodic || _nonbondedMethod == CutoffPeriodic) && r_ij > _taperCutoff) {
RealOpenMM delta = r_ij - _taperCutoff; double delta = r_ij - _taperCutoff;
RealOpenMM taper = 1.0 + delta*delta*delta*(_taperCoefficients[C3] + delta*(_taperCoefficients[C4] + delta*_taperCoefficients[C5])); double taper = 1.0 + delta*delta*delta*(_taperCoefficients[C3] + delta*(_taperCoefficients[C4] + delta*_taperCoefficients[C5]));
RealOpenMM dtaper = delta*delta*(3.0*_taperCoefficients[C3] + delta*(4.0*_taperCoefficients[C4] + delta*5.0*_taperCoefficients[C5])); double dtaper = delta*delta*(3.0*_taperCoefficients[C3] + delta*(4.0*_taperCoefficients[C4] + delta*5.0*_taperCoefficients[C5]));
dEdR = energy*dtaper + dEdR*taper; dEdR = energy*dtaper + dEdR*taper;
energy *= taper; energy *= taper;
} }
dEdR /= r_ij; dEdR /= r_ij;
...@@ -246,16 +230,14 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn(RealOpenMM combinedSigma, R ...@@ -246,16 +230,14 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn(RealOpenMM combinedSigma, R
} }
void AmoebaReferenceVdwForce::setReducedPositions(int numParticles, void AmoebaReferenceVdwForce::setReducedPositions(int numParticles,
const vector<RealVec>& particlePositions, const vector<Vec3>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& reductions, const std::vector<double>& reductions,
std::vector<Vec3>& reducedPositions) const { std::vector<Vec3>& reducedPositions) const {
static const RealOpenMM zero = 0.0;
reducedPositions.resize(numParticles); reducedPositions.resize(numParticles);
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++) { for (unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++) {
if (reductions[ii] != zero) { if (reductions[ii] != 0.0) {
int reductionIndex = indexIVs[ii]; int reductionIndex = indexIVs[ii];
reducedPositions[ii] = Vec3(reductions[ii]*(particlePositions[ii][0] - particlePositions[reductionIndex][0]) + particlePositions[reductionIndex][0], reducedPositions[ii] = Vec3(reductions[ii]*(particlePositions[ii][0] - particlePositions[reductionIndex][0]) + particlePositions[reductionIndex][0],
reductions[ii]*(particlePositions[ii][1] - particlePositions[reductionIndex][1]) + particlePositions[reductionIndex][1], reductions[ii]*(particlePositions[ii][1] - particlePositions[reductionIndex][1]) + particlePositions[reductionIndex][1],
...@@ -266,22 +248,14 @@ void AmoebaReferenceVdwForce::setReducedPositions(int numParticles, ...@@ -266,22 +248,14 @@ void AmoebaReferenceVdwForce::setReducedPositions(int numParticles,
} }
} }
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
const vector<RealVec>& particlePositions, const vector<Vec3>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<double>& sigmas,
const std::vector<RealOpenMM>& epsilons, const std::vector<double>& epsilons,
const std::vector<RealOpenMM>& reductions, const std::vector<double>& reductions,
const std::vector< std::set<int> >& allExclusions, const std::vector< std::set<int> >& allExclusions,
vector<RealVec>& forces) const { vector<Vec3>& forces) const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
// ---------------------------------------------------------------------------------------
// set reduced coordinates // set reduced coordinates
...@@ -297,12 +271,12 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, ...@@ -297,12 +271,12 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
// based on reduction factor // based on reduction factor
// (4) reset exclusion vector // (4) reset exclusion vector
RealOpenMM energy = zero; double energy = 0.0;
std::vector<unsigned int> exclusions(numParticles, 0); std::vector<unsigned int> exclusions(numParticles, 0);
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++) { for (unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++) {
RealOpenMM sigmaI = sigmas[ii]; double sigmaI = sigmas[ii];
RealOpenMM epsilonI = epsilons[ii]; double epsilonI = epsilons[ii];
for (std::set<int>::const_iterator jj = allExclusions[ii].begin(); jj != allExclusions[ii].end(); jj++) { for (std::set<int>::const_iterator jj = allExclusions[ii].begin(); jj != allExclusions[ii].end(); jj++) {
exclusions[*jj] = 1; exclusions[*jj] = 1;
} }
...@@ -310,8 +284,8 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, ...@@ -310,8 +284,8 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
for (unsigned int jj = ii+1; jj < static_cast<unsigned int>(numParticles); jj++) { for (unsigned int jj = ii+1; jj < static_cast<unsigned int>(numParticles); jj++) {
if (exclusions[jj] == 0) { if (exclusions[jj] == 0) {
RealOpenMM combinedSigma = (this->*_combineSigmas)(sigmaI, sigmas[jj]); double combinedSigma = (this->*_combineSigmas)(sigmaI, sigmas[jj]);
RealOpenMM combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj]); double combinedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj]);
Vec3 force; Vec3 force;
energy += calculatePairIxn(combinedSigma, combinedEpsilon, energy += calculatePairIxn(combinedSigma, combinedEpsilon,
...@@ -323,14 +297,14 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, ...@@ -323,14 +297,14 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
forces[ii][1] -= force[1]; forces[ii][1] -= force[1];
forces[ii][2] -= force[2]; forces[ii][2] -= force[2];
} else { } else {
addReducedForce(ii, indexIVs[ii], reductions[ii], -one, force, forces); addReducedForce(ii, indexIVs[ii], reductions[ii], -1.0, force, forces);
} }
if (indexIVs[jj] == jj) { if (indexIVs[jj] == jj) {
forces[jj][0] += force[0]; forces[jj][0] += force[0];
forces[jj][1] += force[1]; forces[jj][1] += force[1];
forces[jj][2] += force[2]; forces[jj][2] += force[2];
} else { } else {
addReducedForce(jj, indexIVs[jj], reductions[jj], one, force, forces); addReducedForce(jj, indexIVs[jj], reductions[jj], 1.0, force, forces);
} }
} }
...@@ -344,22 +318,14 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, ...@@ -344,22 +318,14 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
return energy; return energy;
} }
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
const vector<RealVec>& particlePositions, const vector<Vec3>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<double>& sigmas,
const std::vector<RealOpenMM>& epsilons, const std::vector<double>& epsilons,
const std::vector<RealOpenMM>& reductions, const std::vector<double>& reductions,
const NeighborList& neighborList, const NeighborList& neighborList,
vector<RealVec>& forces) const { vector<Vec3>& forces) const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
// ---------------------------------------------------------------------------------------
// set reduced coordinates // set reduced coordinates
...@@ -372,15 +338,15 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, ...@@ -372,15 +338,15 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
// then call addReducedForce() to apportion force to particle and its covalent partner // then call addReducedForce() to apportion force to particle and its covalent partner
// based on reduction factor // based on reduction factor
RealOpenMM energy = zero; double energy = 0.0;
for (unsigned int ii = 0; ii < neighborList.size(); ii++) { for (unsigned int ii = 0; ii < neighborList.size(); ii++) {
OpenMM::AtomPair pair = neighborList[ii]; OpenMM::AtomPair pair = neighborList[ii];
int siteI = pair.first; int siteI = pair.first;
int siteJ = pair.second; int siteJ = pair.second;
RealOpenMM combinedSigma = (this->*_combineSigmas)(sigmas[siteI], sigmas[siteJ]); double combinedSigma = (this->*_combineSigmas)(sigmas[siteI], sigmas[siteJ]);
RealOpenMM combinedEpsilon = (this->*_combineEpsilons)(epsilons[siteI], epsilons[siteJ]); double combinedEpsilon = (this->*_combineEpsilons)(epsilons[siteI], epsilons[siteJ]);
Vec3 force; Vec3 force;
energy += calculatePairIxn(combinedSigma, combinedEpsilon, energy += calculatePairIxn(combinedSigma, combinedEpsilon,
...@@ -391,14 +357,14 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles, ...@@ -391,14 +357,14 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
forces[siteI][1] -= force[1]; forces[siteI][1] -= force[1];
forces[siteI][2] -= force[2]; forces[siteI][2] -= force[2];
} else { } else {
addReducedForce(siteI, indexIVs[siteI], reductions[siteI], -one, force, forces); addReducedForce(siteI, indexIVs[siteI], reductions[siteI], -1.0, force, forces);
} }
if (indexIVs[siteJ] == siteJ) { if (indexIVs[siteJ] == siteJ) {
forces[siteJ][0] += force[0]; forces[siteJ][0] += force[0];
forces[siteJ][1] += force[1]; forces[siteJ][1] += force[1];
forces[siteJ][2] += force[2]; forces[siteJ][2] += force[2];
} else { } else {
addReducedForce(siteJ, indexIVs[siteJ], reductions[siteJ], one, force, forces); addReducedForce(siteJ, indexIVs[siteJ], reductions[siteJ], 1.0, force, forces);
} }
} }
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.h
View file @ 8469621f
...@@ -25,7 +25,6 @@ ...@@ -25,7 +25,6 @@
#ifndef __AmoebaReferenceVdwForce_H__ #ifndef __AmoebaReferenceVdwForce_H__
#define __AmoebaReferenceVdwForce_H__ #define __AmoebaReferenceVdwForce_H__
#include "RealVec.h"
#include "openmm/Vec3.h" #include "openmm/Vec3.h"
#include "ReferenceNeighborList.h" #include "ReferenceNeighborList.h"
#include <string> #include <string>
...@@ -34,7 +33,7 @@ ...@@ -34,7 +33,7 @@
namespace OpenMM { namespace OpenMM {
class AmoebaReferenceVdwForce; class AmoebaReferenceVdwForce;
typedef RealOpenMM (AmoebaReferenceVdwForce::*CombiningFunction)(RealOpenMM x, RealOpenMM y) const; typedef double (AmoebaReferenceVdwForce::*CombiningFunction)(double x, double y) const;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -178,7 +177,7 @@ public: ...@@ -178,7 +177,7 @@ public:
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void setPeriodicBox(OpenMM::RealVec* vectors); void setPeriodicBox(OpenMM::Vec3* vectors);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -197,12 +196,12 @@ public: ...@@ -197,12 +196,12 @@ public:
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy(int numParticles, const std::vector<OpenMM::RealVec>& particlePositions, double calculateForceAndEnergy(int numParticles, const std::vector<OpenMM::Vec3>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons, const std::vector<double>& sigmas, const std::vector<double>& epsilons,
const std::vector<RealOpenMM>& reductions, const std::vector<double>& reductions,
const std::vector< std::set<int> >& vdwExclusions, const std::vector< std::set<int> >& vdwExclusions,
std::vector<OpenMM::RealVec>& forces) const; std::vector<OpenMM::Vec3>& forces) const;
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -221,12 +220,12 @@ public: ...@@ -221,12 +220,12 @@ public:
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy(int numParticles, const std::vector<OpenMM::RealVec>& particlePositions, double calculateForceAndEnergy(int numParticles, const std::vector<OpenMM::Vec3>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons, const std::vector<double>& sigmas, const std::vector<double>& epsilons,
const std::vector<RealOpenMM>& reductions, const std::vector<double>& reductions,
const NeighborList& neighborList, const NeighborList& neighborList,
std::vector<OpenMM::RealVec>& forces) const; std::vector<OpenMM::Vec3>& forces) const;
private: private:
...@@ -242,18 +241,18 @@ private: ...@@ -242,18 +241,18 @@ private:
double _cutoff; double _cutoff;
double _taperCutoffFactor; double _taperCutoffFactor;
double _taperCutoff; double _taperCutoff;
RealOpenMM _taperCoefficients[3]; double _taperCoefficients[3];
RealVec _periodicBoxVectors[3]; Vec3 _periodicBoxVectors[3];
CombiningFunction _combineSigmas; CombiningFunction _combineSigmas;
RealOpenMM arithmeticSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const; double arithmeticSigmaCombiningRule(double sigmaI, double sigmaJ) const;
RealOpenMM geometricSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const; double geometricSigmaCombiningRule(double sigmaI, double sigmaJ) const;
RealOpenMM cubicMeanSigmaCombiningRule(RealOpenMM sigmaI, RealOpenMM sigmaJ) const; double cubicMeanSigmaCombiningRule(double sigmaI, double sigmaJ) const;
CombiningFunction _combineEpsilons; CombiningFunction _combineEpsilons;
RealOpenMM arithmeticEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const; double arithmeticEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
RealOpenMM geometricEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const; double geometricEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
RealOpenMM harmonicEpsilonCombiningRule(RealOpenMM epsilonI, RealOpenMM epsilonJ) const; double harmonicEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
RealOpenMM hhgEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ) const; double hhgEpsilonCombiningRule(double epsilonI, double epsilonJ) const;
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -272,8 +271,8 @@ private: ...@@ -272,8 +271,8 @@ private:
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void setReducedPositions(int numParticles, const std::vector<RealVec>& particlePositions, void setReducedPositions(int numParticles, const std::vector<Vec3>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<RealOpenMM>& reductions, const std::vector<int>& indexIVs, const std::vector<double>& reductions,
std::vector<Vec3>& reducedPositions) const; std::vector<Vec3>& reducedPositions) const;
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -290,8 +289,8 @@ private: ...@@ -290,8 +289,8 @@ private:
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
void addReducedForce(unsigned int particleI, unsigned int particleIV, void addReducedForce(unsigned int particleI, unsigned int particleIV,
RealOpenMM reduction, RealOpenMM sign, double reduction, double sign,
Vec3& force, std::vector<OpenMM::RealVec>& forces) const; Vec3& force, std::vector<OpenMM::Vec3>& forces) const;
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -317,9 +316,9 @@ private: ...@@ -317,9 +316,9 @@ private:
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
RealOpenMM calculatePairIxn(RealOpenMM combindedSigma, RealOpenMM combindedEpsilon, double calculatePairIxn(double combindedSigma, double combindedEpsilon,
const Vec3& particleIPosition, const Vec3& particleJPosition, const Vec3& particleIPosition, const Vec3& particleJPosition,
Vec3& force) const; Vec3& force) const;
}; };
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceWcaDispersionForce.cpp
View file @ 8469621f
...@@ -22,223 +22,212 @@ ...@@ -22,223 +22,212 @@
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/ */
#ifdef WIN32
#define _USE_MATH_DEFINES // Needed to get M_PI
#endif
#include "AmoebaReferenceForce.h" #include "AmoebaReferenceForce.h"
#include "AmoebaReferenceWcaDispersionForce.h" #include "AmoebaReferenceWcaDispersionForce.h"
#include <cmath>
using std::vector; using std::vector;
using namespace OpenMM; using namespace OpenMM;
AmoebaReferenceWcaDispersionForce::AmoebaReferenceWcaDispersionForce(RealOpenMM epso, RealOpenMM epsh, RealOpenMM rmino, RealOpenMM rminh, AmoebaReferenceWcaDispersionForce::AmoebaReferenceWcaDispersionForce(double epso, double epsh, double rmino, double rminh,
RealOpenMM awater, RealOpenMM shctd, RealOpenMM dispoff, RealOpenMM slevy) : double awater, double shctd, double dispoff, double slevy) :
_epso(epso), _epsh(epsh), _rmino(rmino), _rminh(rminh), _awater(awater), _shctd(shctd), _dispoff(dispoff), _slevy(slevy) { _epso(epso), _epsh(epsh), _rmino(rmino), _rminh(rminh), _awater(awater), _shctd(shctd), _dispoff(dispoff), _slevy(slevy) {
} }
RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn(RealOpenMM radiusI, RealOpenMM radiusK, double AmoebaReferenceWcaDispersionForce::calculatePairIxn(double radiusI, double radiusK,
const RealVec& particleIPosition, const Vec3& particleIPosition,
const RealVec& particleJPosition, const Vec3& particleJPosition,
const RealOpenMM* const intermediateValues, const double* const intermediateValues,
Vec3& force) const { Vec3& force) const {
// --------------------------------------------------------------------------------------- static const double PI = M_PI;
static const RealOpenMM zero = 0.0; double xr = particleIPosition[0] - particleJPosition[0];
static const RealOpenMM one = 1.0; double yr = particleIPosition[1] - particleJPosition[1];
static const RealOpenMM two = 2.0; double zr = particleIPosition[2] - particleJPosition[2];
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 RealOpenMM seven = 7.0;
static const RealOpenMM eight = 8.0;
static const RealOpenMM ten = 10.0;
static const RealOpenMM fortyEight = 48.0;
static const RealOpenMM PI = 3.1415926535897932384;
// --------------------------------------------------------------------------------------- double r2 = xr*xr + yr*yr + zr*zr;
double r = sqrt(r2);
double r3 = r2*r;
RealOpenMM xr = particleIPosition[0] - particleJPosition[0]; double sK = radiusK*_shctd;
RealOpenMM yr = particleIPosition[1] - particleJPosition[1]; double sK2 = sK*sK;
RealOpenMM zr = particleIPosition[2] - particleJPosition[2];
RealOpenMM r2 = xr*xr + yr*yr + zr*zr; double rmixo = intermediateValues[RMIXO];
RealOpenMM r = SQRT(r2); double rmixo7 = intermediateValues[RMIXO7];
RealOpenMM r3 = r2*r;
RealOpenMM sK = radiusK*_shctd; double emixo = intermediateValues[EMIXO];
RealOpenMM sK2 = sK*sK;
RealOpenMM rmixo = intermediateValues[RMIXO]; double rmixh = intermediateValues[RMIXH];
RealOpenMM rmixo7 = intermediateValues[RMIXO7]; double rmixh7 = intermediateValues[RMIXH7];
RealOpenMM emixo = intermediateValues[EMIXO]; double emixh = intermediateValues[EMIXH];
RealOpenMM rmixh = intermediateValues[RMIXH]; double ao = intermediateValues[AO];
RealOpenMM rmixh7 = intermediateValues[RMIXH7]; double ah = intermediateValues[AH];
RealOpenMM emixh = intermediateValues[EMIXH]; double sum = 0.0;
double de = 0.0;
RealOpenMM ao = intermediateValues[AO];
RealOpenMM ah = intermediateValues[AH];
RealOpenMM sum = zero;
RealOpenMM de = zero;
if (radiusI < (r + sK)) { if (radiusI < (r + sK)) {
RealOpenMM rmax = (radiusI > (r - sK)) ? radiusI : (r - sK); double rmax = (radiusI > (r - sK)) ? radiusI : (r - sK);
RealOpenMM lik = rmax; double lik = rmax;
RealOpenMM lik2 = lik*lik; double lik2 = lik*lik;
RealOpenMM lik3 = lik2*lik; double lik3 = lik2*lik;
RealOpenMM lik4 = lik2*lik2; double lik4 = lik2*lik2;
if (lik < rmixo) { if (lik < rmixo) {
RealOpenMM uik = (r + sK) < rmixo ? (r + sK) : rmixo; double uik = (r + sK) < rmixo ? (r + sK) : rmixo;
RealOpenMM uik2 = uik*uik; double uik2 = uik*uik;
RealOpenMM uik3 = uik2*uik; double uik3 = uik2*uik;
RealOpenMM uik4 = uik2*uik2; double uik4 = uik2*uik2;
RealOpenMM term = four*PI/(fortyEight*r)* (three*(lik4-uik4) - eight*r*(lik3-uik3) + six*(r2-sK2)*(lik2-uik2)); double term = 4.0*PI/(48.0*r)* (3.0*(lik4-uik4) - 8.0*r*(lik3-uik3) + 6.0*(r2-sK2)*(lik2-uik2));
RealOpenMM dl; double dl;
if (radiusI > (r - sK)) { if (radiusI > (r - sK)) {
dl = -lik2 + two*(r2 + sK2); dl = -lik2 + 2.0*(r2 + sK2);
dl *= lik2; dl *= lik2;
} else { } else {
dl = -lik3 + four*lik2*r - six*lik*r2 + two*lik*sK2 + four*r*(r2 - sK2); dl = -lik3 + 4.0*lik2*r - 6.0*lik*r2 + 2.0*lik*sK2 + 4.0*r*(r2 - sK2);
dl *= lik; dl *= lik;
} }
RealOpenMM du; double du;
if ((r+sK) > rmixo) { if ((r+sK) > rmixo) {
du = -uik2 + two*(r2 + sK2); du = -uik2 + 2.0*(r2 + sK2);
du *= -uik2; du *= -uik2;
} else { } else {
du = -uik3 + four*uik2*r - six*uik*r2 + two*uik*sK2 + four*r*(r2 - sK2); du = -uik3 + 4.0*uik2*r - 6.0*uik*r2 + 2.0*uik*sK2 + 4.0*r*(r2 - sK2);
du *= -uik; du *= -uik;
} }
de = -emixo*PI*(dl+du)/(four*r2); de = -emixo*PI*(dl+du)/(4.0*r2);
sum += -emixo*term; sum += -emixo*term;
} }
if (lik < rmixh) { if (lik < rmixh) {
RealOpenMM uik = (r + sK) < rmixh ? (r + sK) : rmixh; double uik = (r + sK) < rmixh ? (r + sK) : rmixh;
RealOpenMM uik2 = uik*uik; double uik2 = uik*uik;
RealOpenMM uik3 = uik2*uik; double uik3 = uik2*uik;
RealOpenMM uik4 = uik2*uik2; double uik4 = uik2*uik2;
RealOpenMM term = four*PI / (fortyEight*r)*(three*(lik4-uik4) - eight*r*(lik3-uik3) + six*(r2-sK2)*(lik2-uik2)); double term = 4.0*PI / (48.0*r)*(3.0*(lik4-uik4) - 8.0*r*(lik3-uik3) + 6.0*(r2-sK2)*(lik2-uik2));
RealOpenMM dl; double dl;
if (radiusI > (r-sK)) { if (radiusI > (r-sK)) {
dl = -lik2 + two*(r2 + sK2); dl = -lik2 + 2.0*(r2 + sK2);
dl *= lik2; dl *= lik2;
} else { } else {
dl = -lik3 + four*lik2*r - six*lik*r2 + two*lik*sK2 + four*r*(r2 -sK2); dl = -lik3 + 4.0*lik2*r - 6.0*lik*r2 + 2.0*lik*sK2 + 4.0*r*(r2 -sK2);
dl *= lik; dl *= lik;
} }
RealOpenMM du; double du;
if (r+sK > rmixh) { if (r+sK > rmixh) {
du = -uik2 + two*(r2 + sK2); du = -uik2 + 2.0*(r2 + sK2);
du *= -uik2; du *= -uik2;
} else { } else {
du = -uik3 +four*uik2*r - six*uik*r2 + two*uik*sK2 +four*r*(r2 - sK2); du = -uik3 +4.0*uik2*r - 6.0*uik*r2 + 2.0*uik*sK2 +4.0*r*(r2 - sK2);
du *= -uik; du *= -uik;
} }
de -= two*emixh*PI*(dl+du)/(four*r2); de -= 2.0*emixh*PI*(dl+du)/(4.0*r2);
sum -= two*emixh*term; sum -= 2.0*emixh*term;
} }
RealOpenMM uik = r + sK; double uik = r + sK;
RealOpenMM uik2 = uik * uik; double uik2 = uik * uik;
RealOpenMM uik3 = uik2 * uik; double uik3 = uik2 * uik;
RealOpenMM uik4 = uik2 * uik2; double uik4 = uik2 * uik2;
RealOpenMM uik5 = uik3 * uik2; double uik5 = uik3 * uik2;
RealOpenMM uik6 = uik3 * uik3; double uik6 = uik3 * uik3;
RealOpenMM uik10 = uik5 * uik5; double uik10 = uik5 * uik5;
RealOpenMM uik11 = uik5 * uik6; double uik11 = uik5 * uik6;
RealOpenMM uik12 = uik6 * uik6; double uik12 = uik6 * uik6;
RealOpenMM uik13 = uik10 * uik3; double uik13 = uik10 * uik3;
if (uik > rmixo) { if (uik > rmixo) {
RealOpenMM lik = rmax > rmixo ? rmax : rmixo; double lik = rmax > rmixo ? rmax : rmixo;
RealOpenMM lik2 = lik * lik; double lik2 = lik * lik;
RealOpenMM lik3 = lik2 * lik; double lik3 = lik2 * lik;
RealOpenMM lik4 = lik2 * lik2; double lik4 = lik2 * lik2;
RealOpenMM lik5 = lik2 * lik3; double lik5 = lik2 * lik3;
RealOpenMM lik6 = lik3 * lik3; double lik6 = lik3 * lik3;
RealOpenMM lik10 = lik5 * lik5; double lik10 = lik5 * lik5;
RealOpenMM lik11 = lik5 * lik6; double lik11 = lik5 * lik6;
RealOpenMM lik12 = lik6 * lik6; double lik12 = lik6 * lik6;
RealOpenMM lik13 = lik10 * lik3; double lik13 = lik10 * lik3;
RealOpenMM term = four*PI/(120.0*r*lik5*uik5)*(15.0*uik*lik*r*(uik4-lik4) - ten*uik2*lik2*(uik3-lik3) + six*(sK2-r2)*(uik5-lik5)); double term = 4.0*PI/(120.0*r*lik5*uik5)*(15.0*uik*lik*r*(uik4-lik4) - 10.0*uik2*lik2*(uik3-lik3) + 6.0*(sK2-r2)*(uik5-lik5));
RealOpenMM dl; double dl;
if (radiusI > (r-sK) || rmax < rmixo) { if (radiusI > (r-sK) || rmax < rmixo) {
dl = -five*lik2 + three*(r2 + sK2); dl = -5.0*lik2 + 3.0*(r2 + sK2);
dl /= -lik5; dl /= -lik5;
} else { } else {
dl = five*lik3 - 33.0*lik*r2 - three*lik*sK2 + 15.0*(lik2*r+r3-r*sK2); dl = 5.0*lik3 - 33.0*lik*r2 - 3.0*lik*sK2 + 15.0*(lik2*r+r3-r*sK2);
dl /= lik6; dl /= lik6;
} }
RealOpenMM du = five*uik3 - 33.0*uik*r2 - three*uik*sK2 + 15.0*(uik2*r+r3-r*sK2); double du = 5.0*uik3 - 33.0*uik*r2 - 3.0*uik*sK2 + 15.0*(uik2*r+r3-r*sK2);
du /= -uik6; du /= -uik6;
RealOpenMM idisp = -two*ao*term; double idisp = -2.0*ao*term;
de -= two*ao*PI*(dl + du)/(15.0*r2); de -= 2.0*ao*PI*(dl + du)/(15.0*r2);
term = four*PI/(2640.0*r*lik12*uik12) * (120.0*uik*lik*r*(uik11-lik11) - 66.0*uik2*lik2*(uik10-lik10) + 55.0*(sK2-r2)*(uik12-lik12)); term = 4.0*PI/(2640.0*r*lik12*uik12) * (120.0*uik*lik*r*(uik11-lik11) - 66.0*uik2*lik2*(uik10-lik10) + 55.0*(sK2-r2)*(uik12-lik12));
if (radiusI > (r-sK) || rmax < rmixo) { if (radiusI > (r-sK) || rmax < rmixo) {
dl = -six*lik2 + five*r2 + five*sK2; dl = -6.0*lik2 + 5.0*r2 + 5.0*sK2;
dl /= -lik12; dl /= -lik12;
} else { } else {
dl = six*lik3 - 125.0*lik*r2 - five*lik*sK2 + 60.0*(lik2*r+r3-r*sK2); dl = 6.0*lik3 - 125.0*lik*r2 - 5.0*lik*sK2 + 60.0*(lik2*r+r3-r*sK2);
dl /= lik13; dl /= lik13;
} }
du = six*uik3 - 125.0*uik*r2 - five*uik*sK2 + 60.0*(uik2*r+r3-r*sK2); du = 6.0*uik3 - 125.0*uik*r2 - 5.0*uik*sK2 + 60.0*(uik2*r+r3-r*sK2);
du /= -uik13; du /= -uik13;
de += ao*rmixo7*PI*(dl + du)/(60.0*r2); de += ao*rmixo7*PI*(dl + du)/(60.0*r2);
sum += ao*rmixo7*term + idisp; sum += ao*rmixo7*term + idisp;
} }
if (uik > rmixh) { if (uik > rmixh) {
lik = rmax > rmixh ? rmax : rmixh; lik = rmax > rmixh ? rmax : rmixh;
lik2 = lik * lik; lik2 = lik * lik;
lik3 = lik2 * lik; lik3 = lik2 * lik;
lik4 = lik2 * lik2; lik4 = lik2 * lik2;
RealOpenMM lik5 = lik2 * lik3; double lik5 = lik2 * lik3;
RealOpenMM lik6 = lik3 * lik3; double lik6 = lik3 * lik3;
RealOpenMM lik10 = lik5 * lik5; double lik10 = lik5 * lik5;
RealOpenMM lik11 = lik5 * lik6; double lik11 = lik5 * lik6;
RealOpenMM lik12 = lik6 * lik6; double lik12 = lik6 * lik6;
RealOpenMM lik13 = lik3 * lik10; double lik13 = lik3 * lik10;
RealOpenMM term = four*PI / (120.0*r*lik5*uik5)*(15.0*uik*lik*r*(uik4-lik4) - ten*uik2*lik2*(uik3-lik3) + six*(sK2-r2)*(uik5-lik5)); double term = 4.0*PI / (120.0*r*lik5*uik5)*(15.0*uik*lik*r*(uik4-lik4) - 10.0*uik2*lik2*(uik3-lik3) + 6.0*(sK2-r2)*(uik5-lik5));
RealOpenMM dl; double dl;
if (radiusI > (r-sK) || rmax < rmixh) { if (radiusI > (r-sK) || rmax < rmixh) {
dl = -five*lik2 + three*(r2 + sK2); dl = -5.0*lik2 + 3.0*(r2 + sK2);
dl /= -lik5; dl /= -lik5;
} else { } else {
dl = five*lik3 - 33.0*lik*r2 - three*lik*sK2 + 15.0*(lik2*r+r3-r*sK2); dl = 5.0*lik3 - 33.0*lik*r2 - 3.0*lik*sK2 + 15.0*(lik2*r+r3-r*sK2);
dl /= lik6; dl /= lik6;
} }
RealOpenMM du = five*uik3 - 33.0*uik*r2 - 3.0*uik*sK2 + 15.0*(uik2*r+r3-r*sK2); double du = 5.0*uik3 - 33.0*uik*r2 - 3.0*uik*sK2 + 15.0*(uik2*r+r3-r*sK2);
du /= -uik6; du /= -uik6;
RealOpenMM idisp = -four*ah*term; double idisp = -4.0*ah*term;
de = de - four*ah*PI*(dl + du)/(15.0*r2); de = de - 4.0*ah*PI*(dl + du)/(15.0*r2);
term = four*PI / (2640.0*r*lik12*uik12)* (120.0*uik*lik*r*(uik11-lik11)- 66.0*uik2*lik2*(uik10-lik10)+ 55.0*(sK2-r2)*(uik12-lik12)); term = 4.0*PI / (2640.0*r*lik12*uik12)* (120.0*uik*lik*r*(uik11-lik11)- 66.0*uik2*lik2*(uik10-lik10)+ 55.0*(sK2-r2)*(uik12-lik12));
if (radiusI > (r-sK) || rmax < rmixh) { if (radiusI > (r-sK) || rmax < rmixh) {
dl = -six*lik2 + five*r2 + five*sK2; dl = -6.0*lik2 + 5.0*r2 + 5.0*sK2;
dl = -dl / lik12; dl = -dl / lik12;
} else {; } else {;
dl = six*lik3 - 125.0*lik*r2 - five*lik*sK2 + 60.0*(lik2*r+r3-r*sK2); dl = 6.0*lik3 - 125.0*lik*r2 - 5.0*lik*sK2 + 60.0*(lik2*r+r3-r*sK2);
dl = dl / lik13; dl = dl / lik13;
} }
du = six*uik3 - 125.0*uik*r2 - five*uik*sK2 + 60.0*(uik2*r+r3-r*sK2); du = 6.0*uik3 - 125.0*uik*r2 - 5.0*uik*sK2 + 60.0*(uik2*r+r3-r*sK2);
du /= -uik13; du /= -uik13;
RealOpenMM irep = two*ah*rmixh7*term; double irep = 2.0*ah*rmixh7*term;
de += ah*rmixh7*PI*(dl+du)/(30.0*r2); de += ah*rmixh7*PI*(dl+du)/(30.0*r2);
sum += irep + idisp; sum += irep + idisp;
} }
} }
...@@ -254,70 +243,59 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn(RealOpenMM radius ...@@ -254,70 +243,59 @@ RealOpenMM AmoebaReferenceWcaDispersionForce::calculatePairIxn(RealOpenMM radius
} }
RealOpenMM AmoebaReferenceWcaDispersionForce::calculateForceAndEnergy(int numParticles, double AmoebaReferenceWcaDispersionForce::calculateForceAndEnergy(int numParticles,
const vector<RealVec>& particlePositions, const vector<Vec3>& particlePositions,
const std::vector<RealOpenMM>& radii, const std::vector<double>& radii,
const std::vector<RealOpenMM>& epsilons, const std::vector<double>& epsilons,
RealOpenMM totalMaximumDispersionEnergy, double totalMaximumDispersionEnergy,
vector<RealVec>& forces) const { vector<Vec3>& forces) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceWcaDispersionForce::calculateForceAndEnergy";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM four = 4.0;
// ---------------------------------------------------------------------------------------
// loop over all ixns // loop over all ixns
RealOpenMM energy = zero; double energy = 0.0;
RealOpenMM rmino2 = _rmino*_rmino; double rmino2 = _rmino*_rmino;
RealOpenMM rmino3 = rmino2*_rmino; double rmino3 = rmino2*_rmino;
RealOpenMM rminh2 = _rminh*_rminh; double rminh2 = _rminh*_rminh;
RealOpenMM rminh3 = rminh2*_rminh; double rminh3 = rminh2*_rminh;
RealOpenMM intermediateValues[LastIntermediateValueIndex]; double intermediateValues[LastIntermediateValueIndex];
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++) { for (unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++) {
RealOpenMM epsi = epsilons[ii]; double epsi = epsilons[ii];
RealOpenMM rmini = radii[ii]; double rmini = radii[ii];
RealOpenMM denominator = SQRT(_epso) + SQRT(epsi); double denominator = sqrt(_epso) + sqrt(epsi);
RealOpenMM emixo = four*_epso*epsi/(denominator*denominator); double emixo = 4.0*_epso*epsi/(denominator*denominator);
intermediateValues[EMIXO] = emixo; intermediateValues[EMIXO] = emixo;
RealOpenMM rminI2 = rmini*rmini; double rminI2 = rmini*rmini;
RealOpenMM rminI3 = rminI2*rmini; double rminI3 = rminI2*rmini;
RealOpenMM rmixo = two*(rmino3 + rminI3) / (rmino2 + rminI2); double rmixo = 2.0*(rmino3 + rminI3) / (rmino2 + rminI2);
intermediateValues[RMIXO] = rmixo; intermediateValues[RMIXO] = rmixo;
RealOpenMM rmixo7 = rmixo*rmixo*rmixo; double rmixo7 = rmixo*rmixo*rmixo;
rmixo7 = rmixo7*rmixo7*rmixo; rmixo7 = rmixo7*rmixo7*rmixo;
intermediateValues[RMIXO7] = rmixo7; intermediateValues[RMIXO7] = rmixo7;
intermediateValues[AO] = emixo*rmixo7; intermediateValues[AO] = emixo*rmixo7;
denominator = SQRT(_epsh) + SQRT(epsi); denominator = sqrt(_epsh) + sqrt(epsi);
RealOpenMM emixh = four*_epsh*epsi/ (denominator*denominator); double emixh = 4.0*_epsh*epsi/ (denominator*denominator);
intermediateValues[EMIXH] = emixh; intermediateValues[EMIXH] = emixh;
RealOpenMM rmixh = two * (rminh3 + rminI3) / (rminh2 + rminI2); double rmixh = 2.0 * (rminh3 + rminI3) / (rminh2 + rminI2);
intermediateValues[RMIXH] = rmixh; intermediateValues[RMIXH] = rmixh;
RealOpenMM rmixh7 = rmixh*rmixh*rmixh; double rmixh7 = rmixh*rmixh*rmixh;
rmixh7 = rmixh7*rmixh7*rmixh; rmixh7 = rmixh7*rmixh7*rmixh;
intermediateValues[RMIXH7] = rmixh7; intermediateValues[RMIXH7] = rmixh7;
intermediateValues[AH] = emixh*rmixh7; intermediateValues[AH] = emixh*rmixh7;
for (unsigned int jj = 0; jj < static_cast<unsigned int>(numParticles); jj++) { for (unsigned int jj = 0; jj < static_cast<unsigned int>(numParticles); jj++) {
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceWcaDispersionForce.h
View file @ 8469621f
...@@ -25,7 +25,6 @@ ...@@ -25,7 +25,6 @@
#ifndef __AmoebaReferenceWcaDispersionForce_H__ #ifndef __AmoebaReferenceWcaDispersionForce_H__
#define __AmoebaReferenceWcaDispersionForce_H__ #define __AmoebaReferenceWcaDispersionForce_H__
#include "RealVec.h"
#include "openmm/Vec3.h" #include "openmm/Vec3.h"
#include <string> #include <string>
#include <vector> #include <vector>
...@@ -51,8 +50,8 @@ public: ...@@ -51,8 +50,8 @@ public:
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
AmoebaReferenceWcaDispersionForce(RealOpenMM epso, RealOpenMM epsh, RealOpenMM rmino, RealOpenMM rminh, AmoebaReferenceWcaDispersionForce(double epso, double epsh, double rmino, double rminh,
RealOpenMM awater, RealOpenMM shctd, RealOpenMM dispoff, RealOpenMM slevy); double awater, double shctd, double dispoff, double slevy);
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
...@@ -77,20 +76,20 @@ public: ...@@ -77,20 +76,20 @@ public:
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy(int numParticles, const std::vector<OpenMM::RealVec>& particlePositions, double calculateForceAndEnergy(int numParticles, const std::vector<OpenMM::Vec3>& particlePositions,
const std::vector<RealOpenMM>& radii, const std::vector<double>& radii,
const std::vector<RealOpenMM>& epsilons, const std::vector<double>& epsilons,
RealOpenMM totalMaximumDispersionEnergy, std::vector<OpenMM::RealVec>& forces) const; double totalMaximumDispersionEnergy, std::vector<OpenMM::Vec3>& forces) const;
private: private:
RealOpenMM _epso; double _epso;
RealOpenMM _epsh; double _epsh;
RealOpenMM _rmino; double _rmino;
RealOpenMM _rminh; double _rminh;
RealOpenMM _awater; double _awater;
RealOpenMM _shctd; double _shctd;
RealOpenMM _dispoff; double _dispoff;
RealOpenMM _slevy; double _slevy;
enum { EMIXO, RMIXO, RMIXO7, AO, EMIXH, RMIXH, RMIXH7, AH, LastIntermediateValueIndex }; enum { EMIXO, RMIXO, RMIXO7, AO, EMIXH, RMIXH, RMIXH7, AH, LastIntermediateValueIndex };
...@@ -109,10 +108,10 @@ private: ...@@ -109,10 +108,10 @@ private:
--------------------------------------------------------------------------------------- */ --------------------------------------------------------------------------------------- */
RealOpenMM calculatePairIxn(RealOpenMM radiusI, RealOpenMM radiusJ, double calculatePairIxn(double radiusI, double radiusJ,
const OpenMM::RealVec& particleIPosition, const OpenMM::RealVec& particleJPosition, const OpenMM::Vec3& particleIPosition, const OpenMM::Vec3& particleJPosition,
const RealOpenMM* const intermediateValues, const double* const intermediateValues,
Vec3& force) const; Vec3& force) const;
}; };
......
plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp
View file @ 8469621f
...@@ -41,19 +41,19 @@ ...@@ -41,19 +41,19 @@
using namespace OpenMM; using namespace OpenMM;
using namespace std; using namespace std;
static vector<RealVec>& extractPositions(ContextImpl& context) { static vector<Vec3>& extractPositions(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData()); 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()); 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()); ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<RealVec>*) data->forces); return *((vector<Vec3>*) data->forces);
} }
static ReferenceConstraints& extractConstraints(ContextImpl& context) { static ReferenceConstraints& extractConstraints(ContextImpl& context) {
...@@ -64,13 +64,13 @@ static ReferenceConstraints& extractConstraints(ContextImpl& context) { ...@@ -64,13 +64,13 @@ static ReferenceConstraints& extractConstraints(ContextImpl& context) {
static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& inverseMasses, double timeShift) { static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& inverseMasses, double timeShift) {
const System& system = context.getSystem(); const System& system = context.getSystem();
int numParticles = system.getNumParticles(); int numParticles = system.getNumParticles();
vector<RealVec>& posData = extractPositions(context); vector<Vec3>& posData = extractPositions(context);
vector<RealVec>& velData = extractVelocities(context); vector<Vec3>& velData = extractVelocities(context);
vector<RealVec>& forceData = extractForces(context); vector<Vec3>& forceData = extractForces(context);
// Compute the shifted velocities. // Compute the shifted velocities.
vector<RealVec> shiftedVel(numParticles); vector<Vec3> shiftedVel(numParticles);
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
if (inverseMasses[i] > 0) if (inverseMasses[i] > 0)
shiftedVel[i] = velData[i]+forceData[i]*(timeShift*inverseMasses[i]); shiftedVel[i] = velData[i]+forceData[i]*(timeShift*inverseMasses[i]);
...@@ -119,8 +119,8 @@ void ReferenceCalcDrudeForceKernel::initialize(const System& system, const Drude ...@@ -119,8 +119,8 @@ void ReferenceCalcDrudeForceKernel::initialize(const System& system, const Drude
} }
double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) { double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& pos = extractPositions(context); vector<Vec3>& pos = extractPositions(context);
vector<RealVec>& force = extractForces(context); vector<Vec3>& force = extractForces(context);
int numParticles = particle.size(); int numParticles = particle.size();
double energy = 0; double energy = 0;
...@@ -133,17 +133,17 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include ...@@ -133,17 +133,17 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
int p3 = particle3[i]; int p3 = particle3[i];
int p4 = particle4[i]; int p4 = particle4[i];
RealOpenMM a1 = (p2 == -1 ? 1 : aniso12[i]); double a1 = (p2 == -1 ? 1 : aniso12[i]);
RealOpenMM a2 = (p3 == -1 || p4 == -1 ? 1 : aniso34[i]); double a2 = (p3 == -1 || p4 == -1 ? 1 : aniso34[i]);
RealOpenMM a3 = 3-a1-a2; double a3 = 3-a1-a2;
RealOpenMM k3 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a3); double k3 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a3);
RealOpenMM k1 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a1) - k3; double k1 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a1) - k3;
RealOpenMM k2 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a2) - k3; double k2 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a2) - k3;
// Compute the isotropic force. // Compute the isotropic force.
RealVec delta = pos[p]-pos[p1]; Vec3 delta = pos[p]-pos[p1];
RealOpenMM r2 = delta.dot(delta); double r2 = delta.dot(delta);
energy += 0.5*k3*r2; energy += 0.5*k3*r2;
force[p] -= delta*k3; force[p] -= delta*k3;
force[p1] += delta*k3; force[p1] += delta*k3;
...@@ -151,13 +151,13 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include ...@@ -151,13 +151,13 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
// Compute the first anisotropic force. // Compute the first anisotropic force.
if (p2 != -1) { if (p2 != -1) {
RealVec dir = pos[p1]-pos[p2]; Vec3 dir = pos[p1]-pos[p2];
RealOpenMM invDist = 1.0/sqrt(dir.dot(dir)); double invDist = 1.0/sqrt(dir.dot(dir));
dir *= invDist; dir *= invDist;
RealOpenMM rprime = dir.dot(delta); double rprime = dir.dot(delta);
energy += 0.5*k1*rprime*rprime; energy += 0.5*k1*rprime*rprime;
RealVec f1 = dir*(k1*rprime); Vec3 f1 = dir*(k1*rprime);
RealVec f2 = (delta-dir*rprime)*(k1*rprime*invDist); Vec3 f2 = (delta-dir*rprime)*(k1*rprime*invDist);
force[p] -= f1; force[p] -= f1;
force[p1] += f1-f2; force[p1] += f1-f2;
force[p2] += f2; force[p2] += f2;
...@@ -166,13 +166,13 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include ...@@ -166,13 +166,13 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
// Compute the second anisotropic force. // Compute the second anisotropic force.
if (p3 != -1 && p4 != -1) { if (p3 != -1 && p4 != -1) {
RealVec dir = pos[p3]-pos[p4]; Vec3 dir = pos[p3]-pos[p4];
RealOpenMM invDist = 1.0/sqrt(dir.dot(dir)); double invDist = 1.0/sqrt(dir.dot(dir));
dir *= invDist; dir *= invDist;
RealOpenMM rprime = dir.dot(delta); double rprime = dir.dot(delta);
energy += 0.5*k2*rprime*rprime; energy += 0.5*k2*rprime*rprime;
RealVec f1 = dir*(k2*rprime); Vec3 f1 = dir*(k2*rprime);
RealVec f2 = (delta-dir*rprime)*(k2*rprime*invDist); Vec3 f2 = (delta-dir*rprime)*(k2*rprime*invDist);
force[p] -= f1; force[p] -= f1;
force[p1] += f1; force[p1] += f1;
force[p3] -= f2; force[p3] -= f2;
...@@ -188,18 +188,18 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include ...@@ -188,18 +188,18 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
int dipole2 = pair2[i]; int dipole2 = pair2[i];
int dipole1Particles[] = {particle[dipole1], particle1[dipole1]}; int dipole1Particles[] = {particle[dipole1], particle1[dipole1]};
int dipole2Particles[] = {particle[dipole2], particle1[dipole2]}; int dipole2Particles[] = {particle[dipole2], particle1[dipole2]};
RealOpenMM uscale = pairThole[i]/pow(polarizability[dipole1]*polarizability[dipole2], 1.0/6.0); double uscale = pairThole[i]/pow(polarizability[dipole1]*polarizability[dipole2], 1.0/6.0);
for (int j = 0; j < 2; j++) for (int j = 0; j < 2; j++)
for (int k = 0; k < 2; k++) { for (int k = 0; k < 2; k++) {
int p1 = dipole1Particles[j]; int p1 = dipole1Particles[j];
int p2 = dipole2Particles[k]; int p2 = dipole2Particles[k];
RealOpenMM chargeProduct = charge[dipole1]*charge[dipole2]*(j == k ? 1 : -1); double chargeProduct = charge[dipole1]*charge[dipole2]*(j == k ? 1 : -1);
RealVec delta = pos[p1]-pos[p2]; Vec3 delta = pos[p1]-pos[p2];
RealOpenMM r = sqrt(delta.dot(delta)); double r = sqrt(delta.dot(delta));
RealOpenMM u = r*uscale; double u = r*uscale;
RealOpenMM screening = 1.0 - (1.0+0.5*u)*exp(-u); double screening = 1.0 - (1.0+0.5*u)*exp(-u);
energy += ONE_4PI_EPS0*chargeProduct*screening/r; energy += ONE_4PI_EPS0*chargeProduct*screening/r;
RealVec f = delta*(ONE_4PI_EPS0*chargeProduct/(r*r))*(screening/r-0.5*(1+u)*exp(-u)*uscale); Vec3 f = delta*(ONE_4PI_EPS0*chargeProduct/(r*r))*(screening/r-0.5*(1+u)*exp(-u)*uscale);
force[p1] += f; force[p1] += f;
force[p2] -= f; force[p2] -= f;
} }
...@@ -257,21 +257,21 @@ void ReferenceIntegrateDrudeLangevinStepKernel::initialize(const System& system, ...@@ -257,21 +257,21 @@ void ReferenceIntegrateDrudeLangevinStepKernel::initialize(const System& system,
} }
void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, const DrudeLangevinIntegrator& integrator) { void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
vector<RealVec>& pos = extractPositions(context); vector<Vec3>& pos = extractPositions(context);
vector<RealVec>& vel = extractVelocities(context); vector<Vec3>& vel = extractVelocities(context);
vector<RealVec>& force = extractForces(context); vector<Vec3>& force = extractForces(context);
// Update velocities of ordinary particles. // Update velocities of ordinary particles.
const RealOpenMM vscale = exp(-integrator.getStepSize()*integrator.getFriction()); const double vscale = exp(-integrator.getStepSize()*integrator.getFriction());
const RealOpenMM fscale = (1-vscale)/integrator.getFriction(); const double fscale = (1-vscale)/integrator.getFriction();
const RealOpenMM kT = BOLTZ*integrator.getTemperature(); const double kT = BOLTZ*integrator.getTemperature();
const RealOpenMM noisescale = sqrt(2*kT*integrator.getFriction())*sqrt(0.5*(1-vscale*vscale)/integrator.getFriction()); const double noisescale = sqrt(2*kT*integrator.getFriction())*sqrt(0.5*(1-vscale*vscale)/integrator.getFriction());
for (int i = 0; i < (int) normalParticles.size(); i++) { for (int i = 0; i < (int) normalParticles.size(); i++) {
int index = normalParticles[i]; int index = normalParticles[i];
RealOpenMM invMass = particleInvMass[index]; double invMass = particleInvMass[index];
if (invMass != 0.0) { if (invMass != 0.0) {
RealOpenMM sqrtInvMass = sqrt(invMass); double sqrtInvMass = sqrt(invMass);
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
vel[index][j] = vscale*vel[index][j] + fscale*invMass*force[index][j] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber(); vel[index][j] = vscale*vel[index][j] + fscale*invMass*force[index][j] + noisescale*sqrtInvMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
} }
...@@ -279,21 +279,21 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co ...@@ -279,21 +279,21 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
// Update velocities of Drude particle pairs. // Update velocities of Drude particle pairs.
const RealOpenMM vscaleDrude = exp(-integrator.getStepSize()*integrator.getDrudeFriction()); const double vscaleDrude = exp(-integrator.getStepSize()*integrator.getDrudeFriction());
const RealOpenMM fscaleDrude = (1-vscaleDrude)/integrator.getDrudeFriction(); const double fscaleDrude = (1-vscaleDrude)/integrator.getDrudeFriction();
const RealOpenMM kTDrude = BOLTZ*integrator.getDrudeTemperature(); const double kTDrude = BOLTZ*integrator.getDrudeTemperature();
const RealOpenMM noisescaleDrude = sqrt(2*kTDrude*integrator.getDrudeFriction())*sqrt(0.5*(1-vscaleDrude*vscaleDrude)/integrator.getDrudeFriction()); const double noisescaleDrude = sqrt(2*kTDrude*integrator.getDrudeFriction())*sqrt(0.5*(1-vscaleDrude*vscaleDrude)/integrator.getDrudeFriction());
for (int i = 0; i < (int) pairParticles.size(); i++) { for (int i = 0; i < (int) pairParticles.size(); i++) {
int p1 = pairParticles[i].first; int p1 = pairParticles[i].first;
int p2 = pairParticles[i].second; int p2 = pairParticles[i].second;
RealOpenMM mass1fract = pairInvTotalMass[i]/particleInvMass[p1]; double mass1fract = pairInvTotalMass[i]/particleInvMass[p1];
RealOpenMM mass2fract = pairInvTotalMass[i]/particleInvMass[p2]; double mass2fract = pairInvTotalMass[i]/particleInvMass[p2];
RealOpenMM sqrtInvTotalMass = sqrt(pairInvTotalMass[i]); double sqrtInvTotalMass = sqrt(pairInvTotalMass[i]);
RealOpenMM sqrtInvReducedMass = sqrt(pairInvReducedMass[i]); double sqrtInvReducedMass = sqrt(pairInvReducedMass[i]);
RealVec cmVel = vel[p1]*mass1fract+vel[p2]*mass2fract; Vec3 cmVel = vel[p1]*mass1fract+vel[p2]*mass2fract;
RealVec relVel = vel[p2]-vel[p1]; Vec3 relVel = vel[p2]-vel[p1];
RealVec cmForce = force[p1]+force[p2]; Vec3 cmForce = force[p1]+force[p2];
RealVec relForce = force[p2]*mass1fract - force[p1]*mass2fract; Vec3 relForce = force[p2]*mass1fract - force[p1]*mass2fract;
for (int j = 0; j < 3; j++) { for (int j = 0; j < 3; j++) {
cmVel[j] = vscale*cmVel[j] + fscale*pairInvTotalMass[i]*cmForce[j] + noisescale*sqrtInvTotalMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber(); cmVel[j] = vscale*cmVel[j] + fscale*pairInvTotalMass[i]*cmForce[j] + noisescale*sqrtInvTotalMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
relVel[j] = vscaleDrude*relVel[j] + fscaleDrude*pairInvReducedMass[i]*relForce[j] + noisescaleDrude*sqrtInvReducedMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber(); relVel[j] = vscaleDrude*relVel[j] + fscaleDrude*pairInvReducedMass[i]*relForce[j] + noisescaleDrude*sqrtInvReducedMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
...@@ -305,8 +305,8 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co ...@@ -305,8 +305,8 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
// Update the particle positions. // Update the particle positions.
int numParticles = particleInvMass.size(); int numParticles = particleInvMass.size();
vector<RealVec> xPrime(numParticles); vector<Vec3> xPrime(numParticles);
RealOpenMM dt = integrator.getStepSize(); double dt = integrator.getStepSize();
for (int i = 0; i < numParticles; i++) for (int i = 0; i < numParticles; i++)
if (particleInvMass[i] != 0.0) if (particleInvMass[i] != 0.0)
xPrime[i] = pos[i]+vel[i]*dt; xPrime[i] = pos[i]+vel[i]*dt;
...@@ -317,7 +317,7 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co ...@@ -317,7 +317,7 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
// Record the constrained positions and velocities. // Record the constrained positions and velocities.
RealOpenMM dtInv = 1.0/dt; double dtInv = 1.0/dt;
for (int i = 0; i < numParticles; i++) { for (int i = 0; i < numParticles; i++) {
if (particleInvMass[i] != 0.0) { if (particleInvMass[i] != 0.0) {
vel[i] = (xPrime[i]-pos[i])*dtInv; vel[i] = (xPrime[i]-pos[i])*dtInv;
...@@ -327,31 +327,31 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co ...@@ -327,31 +327,31 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
// Apply hard wall constraints. // Apply hard wall constraints.
const RealOpenMM maxDrudeDistance = integrator.getMaxDrudeDistance(); const double maxDrudeDistance = integrator.getMaxDrudeDistance();
if (maxDrudeDistance > 0) { if (maxDrudeDistance > 0) {
const RealOpenMM hardwallscaleDrude = sqrt(kTDrude); const double hardwallscaleDrude = sqrt(kTDrude);
for (int i = 0; i < (int) pairParticles.size(); i++) { for (int i = 0; i < (int) pairParticles.size(); i++) {
int p1 = pairParticles[i].first; int p1 = pairParticles[i].first;
int p2 = pairParticles[i].second; int p2 = pairParticles[i].second;
RealVec delta = pos[p1]-pos[p2]; Vec3 delta = pos[p1]-pos[p2];
RealOpenMM r = sqrt(delta.dot(delta)); double r = sqrt(delta.dot(delta));
RealOpenMM rInv = 1/r; double rInv = 1/r;
if (rInv*maxDrudeDistance < 1.0) { if (rInv*maxDrudeDistance < 1.0) {
// The constraint has been violated, so make the inter-particle distance "bounce" // The constraint has been violated, so make the inter-particle distance "bounce"
// off the hard wall. // off the hard wall.
if (rInv*maxDrudeDistance < 0.5) if (rInv*maxDrudeDistance < 0.5)
throw OpenMMException("Drude particle moved too far beyond hard wall constraint"); throw OpenMMException("Drude particle moved too far beyond hard wall constraint");
RealVec bondDir = delta*rInv; Vec3 bondDir = delta*rInv;
RealVec vel1 = vel[p1]; Vec3 vel1 = vel[p1];
RealVec vel2 = vel[p2]; Vec3 vel2 = vel[p2];
RealOpenMM mass1 = particleMass[p1]; double mass1 = particleMass[p1];
RealOpenMM mass2 = particleMass[p2]; double mass2 = particleMass[p2];
RealOpenMM deltaR = r-maxDrudeDistance; double deltaR = r-maxDrudeDistance;
RealOpenMM deltaT = dt; double deltaT = dt;
RealOpenMM dotvr1 = vel1.dot(bondDir); double dotvr1 = vel1.dot(bondDir);
RealVec vb1 = bondDir*dotvr1; Vec3 vb1 = bondDir*dotvr1;
RealVec vp1 = vel1-vb1; Vec3 vp1 = vel1-vb1;
if (mass2 == 0) { if (mass2 == 0) {
// The parent particle is massless, so move only the Drude particle. // The parent particle is massless, so move only the Drude particle.
...@@ -360,29 +360,29 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co ...@@ -360,29 +360,29 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
if (deltaT > dt) if (deltaT > dt)
deltaT = dt; deltaT = dt;
dotvr1 = -dotvr1*hardwallscaleDrude/(abs(dotvr1)*sqrt(mass1)); dotvr1 = -dotvr1*hardwallscaleDrude/(abs(dotvr1)*sqrt(mass1));
RealOpenMM dr = -deltaR + deltaT*dotvr1; double dr = -deltaR + deltaT*dotvr1;
pos[p1] += bondDir*dr; pos[p1] += bondDir*dr;
vel[p1] = vp1 + bondDir*dotvr1; vel[p1] = vp1 + bondDir*dotvr1;
} }
else { else {
// Move both particles. // Move both particles.
RealOpenMM invTotalMass = pairInvTotalMass[i]; double invTotalMass = pairInvTotalMass[i];
RealOpenMM dotvr2 = vel2.dot(bondDir); double dotvr2 = vel2.dot(bondDir);
RealVec vb2 = bondDir*dotvr2; Vec3 vb2 = bondDir*dotvr2;
RealVec vp2 = vel2-vb2; Vec3 vp2 = vel2-vb2;
RealOpenMM vbCMass = (mass1*dotvr1 + mass2*dotvr2)*invTotalMass; double vbCMass = (mass1*dotvr1 + mass2*dotvr2)*invTotalMass;
dotvr1 -= vbCMass; dotvr1 -= vbCMass;
dotvr2 -= vbCMass; dotvr2 -= vbCMass;
if (dotvr1 != dotvr2) if (dotvr1 != dotvr2)
deltaT = deltaR/abs(dotvr1-dotvr2); deltaT = deltaR/abs(dotvr1-dotvr2);
if (deltaT > dt) if (deltaT > dt)
deltaT = dt; deltaT = dt;
RealOpenMM vBond = hardwallscaleDrude/sqrt(mass1); double vBond = hardwallscaleDrude/sqrt(mass1);
dotvr1 = -dotvr1*vBond*mass2*invTotalMass/abs(dotvr1); dotvr1 = -dotvr1*vBond*mass2*invTotalMass/abs(dotvr1);
dotvr2 = -dotvr2*vBond*mass1*invTotalMass/abs(dotvr2); dotvr2 = -dotvr2*vBond*mass1*invTotalMass/abs(dotvr2);
RealOpenMM dr1 = -deltaR*mass2*invTotalMass + deltaT*dotvr1; double dr1 = -deltaR*mass2*invTotalMass + deltaT*dotvr1;
RealOpenMM dr2 = deltaR*mass1*invTotalMass + deltaT*dotvr2; double dr2 = deltaR*mass1*invTotalMass + deltaT*dotvr2;
dotvr1 += vbCMass; dotvr1 += vbCMass;
dotvr2 += vbCMass; dotvr2 += vbCMass;
pos[p1] += bondDir*dr1; pos[p1] += bondDir*dr1;
...@@ -419,7 +419,7 @@ void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, cons ...@@ -419,7 +419,7 @@ void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, cons
// Record particle masses. // Record particle masses.
vector<RealOpenMM> particleMass; vector<double> particleMass;
for (int i = 0; i < system.getNumParticles(); i++) { for (int i = 0; i < system.getNumParticles(); i++) {
double mass = system.getParticleMass(i); double mass = system.getParticleMass(i);
particleMass.push_back(mass); particleMass.push_back(mass);
...@@ -433,20 +433,20 @@ void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, cons ...@@ -433,20 +433,20 @@ void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, cons
throw OpenMMException("DrudeSCFIntegrator: Failed to allocate memory"); throw OpenMMException("DrudeSCFIntegrator: Failed to allocate memory");
lbfgs_parameter_init(&minimizerParams); lbfgs_parameter_init(&minimizerParams);
minimizerParams.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE; minimizerParams.linesearch = LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE;
if (sizeof(RealOpenMM) < 8) if (sizeof(double) < 8)
minimizerParams.xtol = 1e-7; minimizerParams.xtol = 1e-7;
} }
void ReferenceIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const DrudeSCFIntegrator& integrator) { void ReferenceIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const DrudeSCFIntegrator& integrator) {
vector<RealVec>& pos = extractPositions(context); vector<Vec3>& pos = extractPositions(context);
vector<RealVec>& vel = extractVelocities(context); vector<Vec3>& vel = extractVelocities(context);
vector<RealVec>& force = extractForces(context); vector<Vec3>& force = extractForces(context);
// Update the positions and velocities. // Update the positions and velocities.
int numParticles = particleInvMass.size(); int numParticles = particleInvMass.size();
vector<RealVec> xPrime(numParticles); vector<Vec3> xPrime(numParticles);
RealOpenMM dt = integrator.getStepSize(); double dt = integrator.getStepSize();
for (int i = 0; i < numParticles; i++) { for (int i = 0; i < numParticles; i++) {
if (particleInvMass[i] != 0.0) { if (particleInvMass[i] != 0.0) {
vel[i] += force[i]*particleInvMass[i]*dt; vel[i] += force[i]*particleInvMass[i]*dt;
...@@ -460,7 +460,7 @@ void ReferenceIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const D ...@@ -460,7 +460,7 @@ void ReferenceIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const D
// Record the constrained positions and velocities. // Record the constrained positions and velocities.
RealOpenMM dtInv = 1.0/dt; double dtInv = 1.0/dt;
for (int i = 0; i < numParticles; i++) { for (int i = 0; i < numParticles; i++) {
if (particleInvMass[i] != 0.0) { if (particleInvMass[i] != 0.0) {
vel[i] = (xPrime[i]-pos[i])*dtInv; vel[i] = (xPrime[i]-pos[i])*dtInv;
...@@ -494,13 +494,13 @@ static lbfgsfloatval_t evaluate(void *instance, const lbfgsfloatval_t *x, lbfgsf ...@@ -494,13 +494,13 @@ static lbfgsfloatval_t evaluate(void *instance, const lbfgsfloatval_t *x, lbfgsf
// Compute the force and energy for this configuration. // Compute the force and energy for this configuration.
vector<RealVec>& pos = extractPositions(context); vector<Vec3>& pos = extractPositions(context);
vector<RealVec>& force = extractForces(context); vector<Vec3>& force = extractForces(context);
for (int i = 0; i < numDrudeParticles; i++) for (int i = 0; i < numDrudeParticles; i++)
pos[drudeParticles[i]] = RealVec(x[3*i], x[3*i+1], x[3*i+2]); pos[drudeParticles[i]] = Vec3(x[3*i], x[3*i+1], x[3*i+2]);
double energy = context.calcForcesAndEnergy(true, true); double energy = context.calcForcesAndEnergy(true, true);
for (int i = 0; i < numDrudeParticles; i++) { for (int i = 0; i < numDrudeParticles; i++) {
RealVec f = force[drudeParticles[i]]; Vec3 f = force[drudeParticles[i]];
g[3*i] = -f[0]; g[3*i] = -f[0];
g[3*i+1] = -f[1]; g[3*i+1] = -f[1];
g[3*i+2] = -f[2]; g[3*i+2] = -f[2];
...@@ -511,11 +511,11 @@ static lbfgsfloatval_t evaluate(void *instance, const lbfgsfloatval_t *x, lbfgsf ...@@ -511,11 +511,11 @@ static lbfgsfloatval_t evaluate(void *instance, const lbfgsfloatval_t *x, lbfgsf
void ReferenceIntegrateDrudeSCFStepKernel::minimize(ContextImpl& context, double tolerance) { void ReferenceIntegrateDrudeSCFStepKernel::minimize(ContextImpl& context, double tolerance) {
// Record the initial positions and determine a normalization constant for scaling the tolerance. // Record the initial positions and determine a normalization constant for scaling the tolerance.
vector<RealVec>& pos = extractPositions(context); vector<Vec3>& pos = extractPositions(context);
int numDrudeParticles = drudeParticles.size(); int numDrudeParticles = drudeParticles.size();
double norm = 0.0; double norm = 0.0;
for (int i = 0; i < numDrudeParticles; i++) { for (int i = 0; i < numDrudeParticles; i++) {
RealVec p = pos[drudeParticles[i]]; Vec3 p = pos[drudeParticles[i]];
minimizerPos[3*i] = p[0]; minimizerPos[3*i] = p[0];
minimizerPos[3*i+1] = p[1]; minimizerPos[3*i+1] = p[1];
minimizerPos[3*i+2] = p[2]; minimizerPos[3*i+2] = p[2];
......
plugins/drude/platforms/reference/src/ReferenceDrudeKernels.h
View file @ 8469621f
...@@ -34,7 +34,7 @@ ...@@ -34,7 +34,7 @@
#include "ReferencePlatform.h" #include "ReferencePlatform.h"
#include "openmm/DrudeKernels.h" #include "openmm/DrudeKernels.h"
#include "RealVec.h" #include "openmm/Vec3.h"
#include "lbfgs.h" #include "lbfgs.h"
#include <utility> #include <utility>
#include <vector> #include <vector>
......
plugins/rpmd/platforms/cuda/tests/TestCudaRpmd.cpp
View file @ 8469621f
...@@ -78,7 +78,7 @@ void testFreeParticles() { ...@@ -78,7 +78,7 @@ void testFreeParticles() {
integ.step(1000); integ.step(1000);
vector<double> ke(numCopies, 0.0); vector<double> ke(numCopies, 0.0);
vector<double> rg(numParticles, 0.0); vector<double> rg(numParticles, 0.0);
const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12); const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
for (int i = 0; i < numSteps; i++) { for (int i = 0; i < numSteps; i++) {
integ.step(1); integ.step(1);
vector<State> state(numCopies); vector<State> state(numCopies);
......
plugins/rpmd/platforms/opencl/tests/TestOpenCLRpmd.cpp
View file @ 8469621f
...@@ -78,7 +78,7 @@ void testFreeParticles() { ...@@ -78,7 +78,7 @@ void testFreeParticles() {
integ.step(1000); integ.step(1000);
vector<double> ke(numCopies, 0.0); vector<double> ke(numCopies, 0.0);
vector<double> rg(numParticles, 0.0); vector<double> rg(numParticles, 0.0);
const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12); const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
for (int i = 0; i < numSteps; i++) { for (int i = 0; i < numSteps; i++) {
integ.step(1); integ.step(1);
vector<State> state(numCopies); vector<State> state(numCopies);
...@@ -171,7 +171,7 @@ void testParaHydrogen() { ...@@ -171,7 +171,7 @@ void testParaHydrogen() {
vector<int> counts(numBins, 0); vector<int> counts(numBins, 0);
const double invBoxSize = 1.0/boxSize; const double invBoxSize = 1.0/boxSize;
double meanKE = 0.0; double meanKE = 0.0;
const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12); const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
for (int step = 0; step < numSteps; step++) { for (int step = 0; step < numSteps; step++) {
integ.step(20); integ.step(20);
vector<State> states(numCopies); vector<State> states(numCopies);
......
plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.cpp
View file @ 8469621f
...@@ -37,19 +37,19 @@ ...@@ -37,19 +37,19 @@
using namespace OpenMM; using namespace OpenMM;
using namespace std; using namespace std;
static vector<RealVec>& extractPositions(ContextImpl& context) { static vector<Vec3>& extractPositions(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData()); 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()); 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()); ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<RealVec>*) data->forces); return *((vector<Vec3>*) data->forces);
} }
ReferenceIntegrateRPMDStepKernel::~ReferenceIntegrateRPMDStepKernel() { ReferenceIntegrateRPMDStepKernel::~ReferenceIntegrateRPMDStepKernel() {
...@@ -113,12 +113,12 @@ void ReferenceIntegrateRPMDStepKernel::initialize(const System& system, const RP ...@@ -113,12 +113,12 @@ void ReferenceIntegrateRPMDStepKernel::initialize(const System& system, const RP
void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDIntegrator& integrator, bool forcesAreValid) { void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDIntegrator& integrator, bool forcesAreValid) {
const int numCopies = positions.size(); const int numCopies = positions.size();
const int numParticles = positions[0].size(); const int numParticles = positions[0].size();
const RealOpenMM dt = integrator.getStepSize(); const double dt = integrator.getStepSize();
const RealOpenMM halfdt = 0.5*dt; const double halfdt = 0.5*dt;
const System& system = context.getSystem(); const System& system = context.getSystem();
vector<RealVec>& pos = extractPositions(context); vector<Vec3>& pos = extractPositions(context);
vector<RealVec>& vel = extractVelocities(context); vector<Vec3>& vel = extractVelocities(context);
vector<RealVec>& f = extractForces(context); vector<Vec3>& f = extractForces(context);
// Loop over copies and compute the force on each one. // Loop over copies and compute the force on each one.
...@@ -129,17 +129,17 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI ...@@ -129,17 +129,17 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
vector<t_complex> v(numCopies); vector<t_complex> v(numCopies);
vector<t_complex> q(numCopies); vector<t_complex> q(numCopies);
const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12); const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
const RealOpenMM scale = 1.0/sqrt((RealOpenMM) numCopies); const double scale = 1.0/sqrt((double) numCopies);
const RealOpenMM nkT = numCopies*BOLTZ*integrator.getTemperature(); const double nkT = numCopies*BOLTZ*integrator.getTemperature();
const RealOpenMM twown = 2.0*nkT/hbar; const double twown = 2.0*nkT/hbar;
const RealOpenMM c1_0 = exp(-halfdt*integrator.getFriction()); const double c1_0 = exp(-halfdt*integrator.getFriction());
const RealOpenMM c2_0 = sqrt(1.0-c1_0*c1_0); const double c2_0 = sqrt(1.0-c1_0*c1_0);
if (integrator.getApplyThermostat()) { if (integrator.getApplyThermostat()) {
for (int particle = 0; particle < numParticles; particle++) { for (int particle = 0; particle < numParticles; particle++) {
if (system.getParticleMass(particle) == 0.0) if (system.getParticleMass(particle) == 0.0)
continue; continue;
const RealOpenMM c3_0 = c2_0*sqrt(nkT/system.getParticleMass(particle)); const double c3_0 = c2_0*sqrt(nkT/system.getParticleMass(particle));
for (int component = 0; component < 3; component++) { for (int component = 0; component < 3; component++) {
for (int k = 0; k < numCopies; k++) for (int k = 0; k < numCopies; k++)
v[k] = t_complex(scale*velocities[k][particle][component], 0.0); v[k] = t_complex(scale*velocities[k][particle][component], 0.0);
...@@ -153,12 +153,12 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI ...@@ -153,12 +153,12 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
for (int k = 1; k <= numCopies/2; k++) { for (int k = 1; k <= numCopies/2; k++) {
const bool isCenter = (numCopies%2 == 0 && k == numCopies/2); const bool isCenter = (numCopies%2 == 0 && k == numCopies/2);
const RealOpenMM wk = twown*sin(k*M_PI/numCopies); const double wk = twown*sin(k*M_PI/numCopies);
const RealOpenMM c1 = exp(-2.0*wk*halfdt); const double c1 = exp(-2.0*wk*halfdt);
const RealOpenMM c2 = sqrt((1.0-c1*c1)/2) * (isCenter ? sqrt(2.0) : 1.0); const double c2 = sqrt((1.0-c1*c1)/2) * (isCenter ? sqrt(2.0) : 1.0);
const RealOpenMM c3 = c2*sqrt(nkT/system.getParticleMass(particle)); const double c3 = c2*sqrt(nkT/system.getParticleMass(particle));
RealOpenMM rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber(); double rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
RealOpenMM rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber()); double rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
v[k] = v[k]*c1 + t_complex(rand1, rand2); v[k] = v[k]*c1 + t_complex(rand1, rand2);
if (k < numCopies-k) if (k < numCopies-k)
v[numCopies-k] = v[numCopies-k]*c1 + t_complex(rand1, -rand2); v[numCopies-k] = v[numCopies-k]*c1 + t_complex(rand1, -rand2);
...@@ -191,11 +191,11 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI ...@@ -191,11 +191,11 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
fftpack_exec_1d(fft, FFTPACK_FORWARD, &v[0], &v[0]); fftpack_exec_1d(fft, FFTPACK_FORWARD, &v[0], &v[0]);
q[0] += v[0]*dt; q[0] += v[0]*dt;
for (int k = 1; k < numCopies; k++) { for (int k = 1; k < numCopies; k++) {
const RealOpenMM wk = twown*sin(k*M_PI/numCopies); const double wk = twown*sin(k*M_PI/numCopies);
const RealOpenMM wt = wk*dt; const double wt = wk*dt;
const RealOpenMM coswt = cos(wt); const double coswt = cos(wt);
const RealOpenMM sinwt = sin(wt); const double sinwt = sin(wt);
const RealOpenMM wm = wk*system.getParticleMass(particle); const double wm = wk*system.getParticleMass(particle);
const t_complex vprime = v[k]*coswt - q[k]*(wk*sinwt); // Advance velocity from t to t+dt const t_complex vprime = v[k]*coswt - q[k]*(wk*sinwt); // Advance velocity from t to t+dt
q[k] = v[k]*(sinwt/wk) + q[k]*coswt; // Advance position from t to t+dt q[k] = v[k]*(sinwt/wk) + q[k]*coswt; // Advance position from t to t+dt
v[k] = vprime; v[k] = vprime;
...@@ -226,7 +226,7 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI ...@@ -226,7 +226,7 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
for (int particle = 0; particle < numParticles; particle++) { for (int particle = 0; particle < numParticles; particle++) {
if (system.getParticleMass(particle) == 0.0) if (system.getParticleMass(particle) == 0.0)
continue; continue;
const RealOpenMM c3_0 = c2_0*sqrt(nkT/system.getParticleMass(particle)); const double c3_0 = c2_0*sqrt(nkT/system.getParticleMass(particle));
for (int component = 0; component < 3; component++) { for (int component = 0; component < 3; component++) {
for (int k = 0; k < numCopies; k++) for (int k = 0; k < numCopies; k++)
v[k] = t_complex(scale*velocities[k][particle][component], 0.0); v[k] = t_complex(scale*velocities[k][particle][component], 0.0);
...@@ -240,12 +240,12 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI ...@@ -240,12 +240,12 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
for (int k = 1; k <= numCopies/2; k++) { for (int k = 1; k <= numCopies/2; k++) {
const bool isCenter = (numCopies%2 == 0 && k == numCopies/2); const bool isCenter = (numCopies%2 == 0 && k == numCopies/2);
const RealOpenMM wk = twown*sin(k*M_PI/numCopies); const double wk = twown*sin(k*M_PI/numCopies);
const RealOpenMM c1 = exp(-2.0*wk*halfdt); const double c1 = exp(-2.0*wk*halfdt);
const RealOpenMM c2 = sqrt((1.0-c1*c1)/2) * (isCenter ? sqrt(2.0) : 1.0); const double c2 = sqrt((1.0-c1*c1)/2) * (isCenter ? sqrt(2.0) : 1.0);
const RealOpenMM c3 = c2*sqrt(nkT/system.getParticleMass(particle)); const double c3 = c2*sqrt(nkT/system.getParticleMass(particle));
RealOpenMM rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber(); double rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
RealOpenMM rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber()); double rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
v[k] = v[k]*c1 + t_complex(rand1, rand2); v[k] = v[k]*c1 + t_complex(rand1, rand2);
if (k < numCopies-k) if (k < numCopies-k)
v[numCopies-k] = v[numCopies-k]*c1 + t_complex(rand1, -rand2); v[numCopies-k] = v[numCopies-k]*c1 + t_complex(rand1, -rand2);
...@@ -265,9 +265,9 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI ...@@ -265,9 +265,9 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const RPMDIntegrator& integrator) { void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const RPMDIntegrator& integrator) {
const int totalCopies = positions.size(); const int totalCopies = positions.size();
const int numParticles = positions[0].size(); const int numParticles = positions[0].size();
vector<RealVec>& pos = extractPositions(context); vector<Vec3>& pos = extractPositions(context);
vector<RealVec>& vel = extractVelocities(context); vector<Vec3>& vel = extractVelocities(context);
vector<RealVec>& f = extractForces(context); vector<Vec3>& f = extractForces(context);
// Compute forces from all groups that didn't have a specified contraction. // Compute forces from all groups that didn't have a specified contraction.
...@@ -298,7 +298,7 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const ...@@ -298,7 +298,7 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const
// Find the contracted positions. // Find the contracted positions.
vector<t_complex> q(totalCopies); vector<t_complex> q(totalCopies);
const RealOpenMM scale1 = 1.0/totalCopies; const double scale1 = 1.0/totalCopies;
for (int particle = 0; particle < numParticles; particle++) { for (int particle = 0; particle < numParticles; particle++) {
for (int component = 0; component < 3; component++) { for (int component = 0; component < 3; component++) {
// Transform to the frequency domain, set high frequency components to zero, and transform back. // Transform to the frequency domain, set high frequency components to zero, and transform back.
...@@ -329,7 +329,7 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const ...@@ -329,7 +329,7 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const
// Apply the forces to the original copies. // Apply the forces to the original copies.
const RealOpenMM scale2 = 1.0/copies; const double scale2 = 1.0/copies;
for (int particle = 0; particle < numParticles; particle++) { for (int particle = 0; particle < numParticles; particle++) {
for (int component = 0; component < 3; component++) { for (int component = 0; component < 3; component++) {
// Transform to the frequency domain, pad with zeros, and transform back. // Transform to the frequency domain, pad with zeros, and transform back.
...@@ -355,12 +355,12 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const ...@@ -355,12 +355,12 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const
double ReferenceIntegrateRPMDStepKernel::computeKineticEnergy(ContextImpl& context, const RPMDIntegrator& integrator) { double ReferenceIntegrateRPMDStepKernel::computeKineticEnergy(ContextImpl& context, const RPMDIntegrator& integrator) {
const System& system = context.getSystem(); const System& system = context.getSystem();
int numParticles = system.getNumParticles(); int numParticles = system.getNumParticles();
vector<RealVec>& velData = extractVelocities(context); vector<Vec3>& velData = extractVelocities(context);
double energy = 0.0; double energy = 0.0;
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
double mass = system.getParticleMass(i); double mass = system.getParticleMass(i);
if (mass > 0) { if (mass > 0) {
RealVec v = velData[i]; Vec3 v = velData[i];
energy += mass*(v.dot(v)); energy += mass*(v.dot(v));
} }
} }
......
plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.h
View file @ 8469621f
...@@ -34,7 +34,7 @@ ...@@ -34,7 +34,7 @@
#include "ReferencePlatform.h" #include "ReferencePlatform.h"
#include "openmm/RpmdKernels.h" #include "openmm/RpmdKernels.h"
#include "RealVec.h" #include "openmm/Vec3.h"
#include "fftpack.h" #include "fftpack.h"
namespace OpenMM { namespace OpenMM {
...@@ -86,11 +86,11 @@ public: ...@@ -86,11 +86,11 @@ public:
void copyToContext(int copy, ContextImpl& context); void copyToContext(int copy, ContextImpl& context);
private: private:
void computeForces(ContextImpl& context, const RPMDIntegrator& integrator); void computeForces(ContextImpl& context, const RPMDIntegrator& integrator);
std::vector<std::vector<RealVec> > positions; std::vector<std::vector<Vec3> > positions;
std::vector<std::vector<RealVec> > velocities; std::vector<std::vector<Vec3> > velocities;
std::vector<std::vector<RealVec> > forces; std::vector<std::vector<Vec3> > forces;
std::vector<std::vector<RealVec> > contractedPositions; std::vector<std::vector<Vec3> > contractedPositions;
std::vector<std::vector<RealVec> > contractedForces; std::vector<std::vector<Vec3> > contractedForces;
std::map<int, int> groupsByCopies; std::map<int, int> groupsByCopies;
int groupsNotContracted; int groupsNotContracted;
fftpack* fft; fftpack* fft;
......
plugins/rpmd/platforms/reference/tests/TestReferenceRpmd.cpp
View file @ 8469621f
...@@ -77,7 +77,7 @@ void testFreeParticles() { ...@@ -77,7 +77,7 @@ void testFreeParticles() {
integ.step(1000); integ.step(1000);
vector<double> ke(numCopies, 0.0); vector<double> ke(numCopies, 0.0);
vector<double> rg(numParticles, 0.0); vector<double> rg(numParticles, 0.0);
const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12); const double hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
for (int i = 0; i < numSteps; i++) { for (int i = 0; i < numSteps; i++) {
integ.step(1); integ.step(1);
vector<State> state(numCopies); vector<State> state(numCopies);
......
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