Skip to content

GitLab

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

Reference platform no longer premultiplies charges by sqrt(1/4*pi*eps0).

parent a170606b
Show whitespace changes
Inline Side-by-side
Showing with 17 additions and 164 deletions
platforms/reference/src/ReferenceKernels.cpp
View file @ cf7d9866
......@@ -498,13 +498,12 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
bonded14IndexArray = allocateIntArray(num14, 2);
bonded14ParamArray = allocateRealArray(num14, 3);
particleParamArray = allocateRealArray(numParticles, 3);
RealOpenMM sqrtEps = static_cast<RealOpenMM>( std::sqrt(138.935485) );
for (int i = 0; i < numParticles; ++i) {
double charge, radius, depth;
force.getParticleParameters(i, charge, radius, depth);
particleParamArray[i][0] = static_cast<RealOpenMM>(0.5*radius);
particleParamArray[i][1] = static_cast<RealOpenMM>(2.0*sqrt(depth));
particleParamArray[i][2] = static_cast<RealOpenMM>(charge*sqrtEps);
particleParamArray[i][2] = static_cast<RealOpenMM>(charge);
}
this->exclusions = exclusions;
exclusionArray = new int*[numParticles];
......@@ -523,7 +522,7 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
bonded14IndexArray[i][1] = particle2;
bonded14ParamArray[i][0] = static_cast<RealOpenMM>(radius);
bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge*sqrtEps*sqrtEps);
bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge);
}
nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
......
platforms/reference/src/SimTKReference/PME.cpp
View file @ cf7d9866
......@@ -37,10 +37,6 @@
#include "fftpack.h"
//#define ONE_4PI_EPS0 (33.20636930*4.184) /* Units of kJ/mol and nm */
// In OpenMM, atom charges are already scaled by sqrt(ONE_4PI_EPS0), so we don't need this
#define ONE_4PI_EPS0 (1.0) /* Units of kJ/mol and nm */
typedef int ivec[3];
......
platforms/reference/src/SimTKReference/ReferenceLJCoulomb14.cpp
View file @ cf7d9866
......@@ -84,51 +84,6 @@ ReferenceLJCoulomb14::~ReferenceLJCoulomb14( ){
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate parameters for LJ 1-4 ixn
@param c6 c6
@param c12 c12
@param q1 q1 charge atom 1
@param q2 q2 charge atom 2
@param epsfac epsfac ????????????
@param parameters output parameters:
parameter[0]= c6*c6/c12
parameter[1]= (c12/c6)**1/6
parameter[2]= epsfactor*q1*q2
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulomb14::getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM q2, RealOpenMM epsfac,
RealOpenMM* parameters ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceLJCoulomb14::getDerivedParameters";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM six = 6.0;
static const RealOpenMM oneSixth = one/six;
// ---------------------------------------------------------------------------------------
if( c12 <= zero ){
parameters[0] = one;
parameters[1] = zero;
} else {
parameters[0] = (c6*c6)/c12;
parameters[1] = POW( (c12/c6), oneSixth );
}
parameters[2] = epsfac*q1*q2;
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate LJ 1-4 ixn
......@@ -202,9 +157,9 @@ int ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomC
RealOpenMM dEdR = parameters[1]*( twelve*sig6 - six )*sig6;
if (cutoff)
dEdR += parameters[2]*(inverseR-2.0f*krf*r2);
dEdR += ONE_4PI_EPS0*parameters[2]*(inverseR-2.0f*krf*r2);
else
dEdR += parameters[2]*inverseR;
dEdR += ONE_4PI_EPS0*parameters[2]*inverseR;
dEdR *= inverseR*inverseR;
// accumulate forces
......@@ -217,9 +172,9 @@ int ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomC
RealOpenMM energy = parameters[1]*( sig6 - one )*sig6;
if (cutoff)
energy += parameters[2]*(inverseR+krf*r2-crf);
energy += ONE_4PI_EPS0*parameters[2]*(inverseR+krf*r2-crf);
else
energy += parameters[2]*inverseR;
energy += ONE_4PI_EPS0*parameters[2]*inverseR;
// accumulate energies
......
platforms/reference/src/SimTKReference/ReferenceLJCoulomb14.h
View file @ cf7d9866
......@@ -68,28 +68,6 @@ class ReferenceLJCoulomb14 : public ReferenceBondIxn {
int setUseCutoff( RealOpenMM distance, RealOpenMM solventDielectric );
/**---------------------------------------------------------------------------------------
Calculate parameters for LJ 1-4 ixn
@param c6 c6
@param c12 c12
@param q1 q1 charge atom 1
@param q2 q2 charge atom 2
@param epsfac epsfac ????????????
@param parameters output parameters:
parameter[0]= c6*c6/c12
parameter[1]= (c12/c6)**1/6
parameter[2]= epsfactor*q1*q2
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM q2, RealOpenMM epsfac,
RealOpenMM* parameters ) const;
/**---------------------------------------------------------------------------------------
Calculate Ryckaert-Bellemans bond ixn
......
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
View file @ cf7d9866
......@@ -156,58 +156,6 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){
pme = true;
}
/**---------------------------------------------------------------------------------------
Calculate parameters for LJ Coulomb ixn
@param c6 c6
@param c12 c12
@param q1 q1 charge atom 1
@param epsfac epsfacSqrt ????????????
@param parameters output parameters:
parameter[SigIndex] = 0.5*( (c12/c6)**1/6 ) (sigma/2)
parameter[EpsIndex] = sqrt(c6*c6/c12) (2*sqrt(epsilon))
parameter[QIndex] = epsfactorSqrt*q1
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulombIxn::getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM epsfacSqrt,
RealOpenMM* parameters ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceLJCoulombIxn::getDerivedParameters";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM six = 6.0;
static const RealOpenMM half = 0.5;
static const RealOpenMM oneSixth = one/six;
static const RealOpenMM oneTweleth = half*oneSixth;
// ---------------------------------------------------------------------------------------
if( c12 <= 0.0 ){
parameters[EpsIndex] = zero;
parameters[SigIndex] = half;
} else {
parameters[EpsIndex] = c6*SQRT( one/c12 );
parameters[SigIndex] = POW( (c12/c6), oneSixth );
parameters[SigIndex] *= half;
}
parameters[QIndex] = epsfacSqrt*q1;
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate Ewald ixn
......@@ -245,7 +193,7 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
RealOpenMM factorEwald = -1 / (4*alphaEwald*alphaEwald);
RealOpenMM SQRT_PI = sqrt(PI);
RealOpenMM TWO_PI = 2.0 * PI;
RealOpenMM recipCoeff = (RealOpenMM)(4*PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
RealOpenMM recipCoeff = (RealOpenMM)(ONE_4PI_EPS0*4*PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
RealOpenMM totalSelfEwaldEnergy = 0.0;
RealOpenMM realSpaceEwaldEnergy = 0.0;
......@@ -258,7 +206,7 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
// **************************************************************************************
for( int atomID = 0; atomID < numberOfAtoms; atomID++ ){
RealOpenMM selfEwaldEnergy = atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI;
RealOpenMM selfEwaldEnergy = ONE_4PI_EPS0*atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI;
totalSelfEwaldEnergy -= selfEwaldEnergy;
if( energyByAtom ){
......@@ -421,12 +369,11 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r;
RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
RealOpenMM dEdR = ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
dEdR = (RealOpenMM)(dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
RealOpenMM sig = atomParameters[ii][SigIndex] + atomParameters[jj][SigIndex];
......@@ -446,7 +393,7 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
// accumulate energies
realSpaceEwaldEnergy = (RealOpenMM) (atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR));
realSpaceEwaldEnergy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR));
vdwEnergy = eps*(sig6-one)*sig6;
totalVdwEnergy += vdwEnergy;
......@@ -476,7 +423,7 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r;
RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
RealOpenMM dEdR = ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
dEdR = (RealOpenMM)(dEdR * (erf(alphaR) - 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
// accumulate forces
......@@ -489,7 +436,7 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
// accumulate energies
realSpaceEwaldEnergy = (RealOpenMM) (atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erf(alphaR));
realSpaceEwaldEnergy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erf(alphaR));
totalExclusionEnergy += realSpaceEwaldEnergy;
if( energyByAtom ){
......@@ -635,9 +582,9 @@ int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoo
RealOpenMM eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
RealOpenMM dEdR = eps*( twelve*sig6 - six )*sig6;
if (cutoff)
dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
dEdR += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
else
dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
dEdR += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
dEdR *= inverseR*inverseR;
// accumulate forces
......@@ -654,9 +601,9 @@ int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoo
if( totalEnergy || energyByAtom ) {
if (cutoff)
energy = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf);
energy = ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf);
else
energy = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
energy = ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
energy += eps*(sig6-one)*sig6;
if( totalEnergy )
*totalEnergy += energy;
......
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.h
View file @ cf7d9866
......@@ -144,28 +144,6 @@ class ReferenceLJCoulombIxn : public ReferencePairIxn {
void setUsePME(RealOpenMM alpha, int meshSize[3]);
/**---------------------------------------------------------------------------------------
Calculate parameters for LJ 1-4 ixn
@param c6 c6
@param c12 c12
@param q1 q1 charge atom
@param epsfacSqrt epsfacSqrt (what is this?)
@param parameters output parameters:
parameter[SigIndex] = sqrt(c6*c6/c12)
parameter[EpsIndex] = 0.5*( (c12/c6)**1/6 )
parameter[QIndex] = epsfactorSqrt*q1
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM epsfacSqrt,
RealOpenMM* parameters ) const;
/**---------------------------------------------------------------------------------------
Calculate LJ Coulomb pair ixn
......
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