Skip to content

GitLab

Commit 042ab1aa authored by Peter Eastman's avatar Peter Eastman
Browse files

Use a consistent (and more accurate) value for electric constant everywhere. ... 

Use a consistent (and more accurate) value for electric constant everywhere.  Also cleaned up handling of constants and units in various places.
parent 2f921936
Hide whitespace changes
Inline Side-by-side
Showing with 90 additions and 109 deletions
openmmapi/include/openmm/CustomGBForce.h
View file @ 042ab1aa
...@@ -97,9 +97,9 @@ namespace OpenMM { ...@@ -97,9 +97,9 @@ namespace OpenMM {
* "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePair); * "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePair);
* custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);" * custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
* "psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle); * "psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle);
* custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", * custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B",
* CustomGBForce::SingleParticle); * CustomGBForce::SingleParticle);
* custom->addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;" * custom->addEnergyTerm("-138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
* "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePair); * "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePair);
* </pre></tt> * </pre></tt>
* *
...@@ -353,7 +353,7 @@ public: ...@@ -353,7 +353,7 @@ public:
* value is being calculated, and "2" to indicate the particle it is interacting with. * value is being calculated, and "2" to indicate the particle it is interacting with.
* @param type the method to use for computing this value * @param type the method to use for computing this value
*/ */
void getComputedValueParameters(int index, std::string& name, std::string& expression, ComputationType type) const; void getComputedValueParameters(int index, std::string& name, std::string& expression, ComputationType& type) const;
/** /**
* Set the properties of a computed value. * Set the properties of a computed value.
* *
...@@ -404,7 +404,7 @@ public: ...@@ -404,7 +404,7 @@ public:
* in the pair and "2" to indicate the second particle in the pair. * in the pair and "2" to indicate the second particle in the pair.
* @param type the method to use for computing this value * @param type the method to use for computing this value
*/ */
void getEnergyTermParameters(int index, std::string& expression, ComputationType type) const; void getEnergyTermParameters(int index, std::string& expression, ComputationType& type) const;
/** /**
* Set the properties of a term to the energy computation. * Set the properties of a term to the energy computation.
* *
......
openmmapi/src/CustomGBForce.cpp
View file @ 042ab1aa
...@@ -117,7 +117,7 @@ int CustomGBForce::addComputedValue(const std::string& name, const std::string& ...@@ -117,7 +117,7 @@ int CustomGBForce::addComputedValue(const std::string& name, const std::string&
return computedValues.size()-1; return computedValues.size()-1;
} }
void CustomGBForce::getComputedValueParameters(int index, std::string& name, std::string& expression, ComputationType type) const { void CustomGBForce::getComputedValueParameters(int index, std::string& name, std::string& expression, ComputationType& type) const {
name = computedValues[index].name; name = computedValues[index].name;
expression = computedValues[index].expression; expression = computedValues[index].expression;
type = computedValues[index].type; type = computedValues[index].type;
...@@ -134,7 +134,7 @@ int CustomGBForce::addEnergyTerm(const std::string& expression, ComputationType ...@@ -134,7 +134,7 @@ int CustomGBForce::addEnergyTerm(const std::string& expression, ComputationType
return energyTerms.size()-1; return energyTerms.size()-1;
} }
void CustomGBForce::getEnergyTermParameters(int index, std::string& expression, ComputationType type) const { void CustomGBForce::getEnergyTermParameters(int index, std::string& expression, ComputationType& type) const {
expression = energyTerms[index].expression; expression = energyTerms[index].expression;
type = energyTerms[index].type; type = energyTerms[index].type;
} }
......
platforms/cuda/src/CudaKernels.cpp
View file @ 042ab1aa
...@@ -31,6 +31,7 @@ ...@@ -31,6 +31,7 @@
#include "openmm/internal/NonbondedForceImpl.h" #include "openmm/internal/NonbondedForceImpl.h"
#include "kernels/gputypes.h" #include "kernels/gputypes.h"
#include "kernels/cudaKernels.h" #include "kernels/cudaKernels.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <cmath> #include <cmath>
extern "C" int gpuSetConstants( gpuContext gpu ); extern "C" int gpuSetConstants( gpuContext gpu );
...@@ -414,7 +415,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon ...@@ -414,7 +415,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
} }
} }
data.nonbondedMethod = method; data.nonbondedMethod = method;
gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method); gpuSetCoulombParameters(gpu, (float) ONE_4PI_EPS0, particle, c6, c12, q, symbol, exclusionList, method);
// Compute the Ewald self energy. // Compute the Ewald self energy.
...@@ -444,7 +445,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon ...@@ -444,7 +445,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
q1[i] = (float) charge; q1[i] = (float) charge;
q2[i] = 1.0f; q2[i] = 1.0f;
} }
gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, particle1, particle2, c6, c12, q1, q2); gpuSetLJ14Parameters(gpu, (float) ONE_4PI_EPS0, 1.0f, particle1, particle2, c6, c12, q1, q2);
} }
} }
...@@ -654,8 +655,8 @@ static void initializeIntegration(const System& system, CudaPlatform::PlatformDa ...@@ -654,8 +655,8 @@ static void initializeIntegration(const System& system, CudaPlatform::PlatformDa
gpuSetRbDihedralParameters(gpu, vector<int>(), vector<int>(), vector<int>(), vector<int>(), vector<float>(), vector<float>(), gpuSetRbDihedralParameters(gpu, vector<int>(), vector<int>(), vector<int>(), vector<int>(), vector<float>(), vector<float>(),
vector<float>(), vector<float>(), vector<float>(), vector<float>()); vector<float>(), vector<float>(), vector<float>(), vector<float>());
if (!data.hasNonbonded) { if (!data.hasNonbonded) {
gpuSetCoulombParameters(gpu, 138.935485f, vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<char>(), vector<vector<int> >(), NO_CUTOFF); gpuSetCoulombParameters(gpu, (float) ONE_4PI_EPS0, vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<char>(), vector<vector<int> >(), NO_CUTOFF);
gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, vector<int>(), vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<float>()); gpuSetLJ14Parameters(gpu, (float) ONE_4PI_EPS0, 1.0f, vector<int>(), vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<float>());
} }
// Set masses. // Set masses.
......
platforms/cuda/src/kernels/gpu.cpp
View file @ 042ab1aa
...@@ -122,14 +122,14 @@ static const float forceConversionFactor = 0.4184f; ...@@ -122,14 +122,14 @@ static const float forceConversionFactor = 0.4184f;
//static const float surfaceAreaFactor = -6.0f * PI * 4.0f * 0.0049f * 1000.0f; //static const float surfaceAreaFactor = -6.0f * PI * 4.0f * 0.0049f * 1000.0f;
static const float surfaceAreaFactor = -6.0f*PI*0.0216f*1000.0f*0.4184f; static const float surfaceAreaFactor = -6.0f*PI*0.0216f*1000.0f*0.4184f;
//static const float surfaceAreaFactor = -1.7035573959e+001; //static const float surfaceAreaFactor = -1.7035573959e+001;
//static const float surfaceAreaFactor = -166.02691f; //static const float surfaceAreaFactor = -166.03185f;
//static const float surfaceAreaFactor = 1.0f; //static const float surfaceAreaFactor = 1.0f;
static const float alphaOBC = 1.0f; static const float alphaOBC = 1.0f;
static const float betaOBC = 0.8f; static const float betaOBC = 0.8f;
static const float gammaOBC = 4.85f; static const float gammaOBC = 4.85f;
static const float kcalMolTokJNM = -0.4184f; static const float kcalMolTokJNM = -0.4184f;
static const float electricConstant = -166.02691f; static const float electricConstant = -166.03185f;
static const float defaultInnerDielectric = 1.0f; static const float defaultInnerDielectric = 1.0f;
static const float defaultSolventDielectric = 78.3f; static const float defaultSolventDielectric = 78.3f;
static const float KILO = 1e3; // Thousand static const float KILO = 1e3; // Thousand
......
platforms/cuda/tests/TestCudaCustomNonbondedForce.cpp
View file @ 042ab1aa
...@@ -250,7 +250,7 @@ void testCoulombLennardJones() { ...@@ -250,7 +250,7 @@ void testCoulombLennardJones() {
customSystem.addParticle(1.0); customSystem.addParticle(1.0);
} }
NonbondedForce* standardNonbonded = new NonbondedForce(); NonbondedForce* standardNonbonded = new NonbondedForce();
CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935485*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)"); CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935456*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
customNonbonded->addPerParticleParameter("q"); customNonbonded->addPerParticleParameter("q");
customNonbonded->addPerParticleParameter("sigma"); customNonbonded->addPerParticleParameter("sigma");
customNonbonded->addPerParticleParameter("eps"); customNonbonded->addPerParticleParameter("eps");
......
platforms/cuda/tests/TestCudaNonbondedForce.cpp
View file @ 042ab1aa
...@@ -71,10 +71,10 @@ void testCoulomb() { ...@@ -71,10 +71,10 @@ void testCoulomb() {
context.setPositions(positions); context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces(); const vector<Vec3>& forces = state.getForces();
double force = 138.935485*(-0.75)/4.0; double force = ONE_4PI_EPS0*(-0.75)/4.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL); ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL); ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(138.935485*(-0.75)/2.0, state.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(-0.75)/2.0, state.getPotentialEnergy(), TOL);
} }
void testLJ() { void testLJ() {
...@@ -173,8 +173,8 @@ void testExclusionsAnd14() { ...@@ -173,8 +173,8 @@ void testExclusionsAnd14() {
context2.setPositions(positions); context2.setPositions(positions);
state = context2.getState(State::Forces | State::Energy); state = context2.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces(); const vector<Vec3>& forces2 = state.getForces();
force = 138.935485*4/(r*r); force = ONE_4PI_EPS0*4/(r*r);
energy = 138.935485*4/r; energy = ONE_4PI_EPS0*4/r;
if (i == 3) { if (i == 3) {
force /= 1.2; force /= 1.2;
energy /= 1.2; energy /= 1.2;
...@@ -216,13 +216,13 @@ void testCutoff() { ...@@ -216,13 +216,13 @@ void testCutoff() {
const vector<Vec3>& forces = state.getForces(); const vector<Vec3>& forces = state.getForces();
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0); const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0); const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force1 = 138.935485*(1.0)*(0.25-2.0*krf*2.0); const double force1 = ONE_4PI_EPS0*(1.0)*(0.25-2.0*krf*2.0);
const double force2 = 138.935485*(1.0)*(1.0-2.0*krf*1.0); const double force2 = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL); ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL); ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL); ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL);
const double energy1 = 138.935485*(1.0)*(0.5+krf*4.0-crf); const double energy1 = ONE_4PI_EPS0*(1.0)*(0.5+krf*4.0-crf);
const double energy2 = 138.935485*(1.0)*(1.0+krf*1.0-crf); const double energy2 = ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf);
ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL);
} }
...@@ -308,8 +308,8 @@ void testCutoff14() { ...@@ -308,8 +308,8 @@ void testCutoff14() {
const vector<Vec3>& forces2 = state.getForces(); const vector<Vec3>& forces2 = state.getForces();
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0); const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0); const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
force = 138.935485*q*q*(1.0/(r*r)-2.0*krf*r); force = ONE_4PI_EPS0*q*q*(1.0/(r*r)-2.0*krf*r);
energy = 138.935485*q*q*(1.0/r+krf*r*r-crf); energy = ONE_4PI_EPS0*q*q*(1.0/r+krf*r*r-crf);
if (i == 3) { if (i == 3) {
force /= 1.2; force /= 1.2;
energy /= 1.2; energy /= 1.2;
...@@ -352,11 +352,11 @@ void testPeriodic() { ...@@ -352,11 +352,11 @@ void testPeriodic() {
const double eps = 78.3; const double eps = 78.3;
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0); const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0); const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force = 138.935485*(1.0)*(1.0-2.0*krf*1.0); const double force = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL); ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL); ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL); ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
} }
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ 042ab1aa
...@@ -35,17 +35,12 @@ ...@@ -35,17 +35,12 @@
#include "OpenCLNonbondedUtilities.h" #include "OpenCLNonbondedUtilities.h"
#include "lepton/Parser.h" #include "lepton/Parser.h"
#include "lepton/ParsedExpression.h" #include "lepton/ParsedExpression.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <cmath> #include <cmath>
using namespace OpenMM; using namespace OpenMM;
using namespace std; using namespace std;
static const double KILO = 1e3; // Thousand
static const double BOLTZMANN = 1.380658e-23; // (J/K)
static const double AVOGADRO = 6.0221367e23; // ()
static const double RGAS = BOLTZMANN*AVOGADRO; // (J/(mol K))
static const double BOLTZ = (RGAS/KILO); // (kJ/(mol K))
static string doubleToString(double value) { static string doubleToString(double value) {
stringstream s; stringstream s;
s.precision(8); s.precision(8);
...@@ -764,8 +759,8 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -764,8 +759,8 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
defines["EWALD_ALPHA"] = doubleToString(alpha); defines["EWALD_ALPHA"] = doubleToString(alpha);
defines["TWO_OVER_SQRT_PI"] = doubleToString(2.0/sqrt(M_PI)); defines["TWO_OVER_SQRT_PI"] = doubleToString(2.0/sqrt(M_PI));
defines["USE_EWALD"] = "1"; defines["USE_EWALD"] = "1";
double selfEnergyScale = 138.935485*alpha/std::sqrt(M_PI); double selfEnergyScale = ONE_4PI_EPS0*alpha/std::sqrt(M_PI);
ewaldSelfEnergy = - 138.935485*alpha*sumSquaredCharges/std::sqrt(M_PI); ewaldSelfEnergy = -ONE_4PI_EPS0*alpha*sumSquaredCharges/std::sqrt(M_PI);
// Create the reciprocal space kernels. // Create the reciprocal space kernels.
...@@ -777,7 +772,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -777,7 +772,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
replacements["RECIPROCAL_BOX_SIZE_X"] = doubleToString(2.0*M_PI/boxVectors[0][0]); replacements["RECIPROCAL_BOX_SIZE_X"] = doubleToString(2.0*M_PI/boxVectors[0][0]);
replacements["RECIPROCAL_BOX_SIZE_Y"] = doubleToString(2.0*M_PI/boxVectors[1][1]); replacements["RECIPROCAL_BOX_SIZE_Y"] = doubleToString(2.0*M_PI/boxVectors[1][1]);
replacements["RECIPROCAL_BOX_SIZE_Z"] = doubleToString(2.0*M_PI/boxVectors[2][2]); replacements["RECIPROCAL_BOX_SIZE_Z"] = doubleToString(2.0*M_PI/boxVectors[2][2]);
replacements["RECIPROCAL_COEFFICIENT"] = doubleToString(138.935485*4*M_PI/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2])); replacements["RECIPROCAL_COEFFICIENT"] = doubleToString(ONE_4PI_EPS0*4*M_PI/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]));
replacements["EXP_COEFFICIENT"] = doubleToString(-1.0/(4.0*alpha*alpha)); replacements["EXP_COEFFICIENT"] = doubleToString(-1.0/(4.0*alpha*alpha));
cl::Program program = cl.createProgram(cl.loadSourceFromFile("ewald.cl"), replacements); cl::Program program = cl.createProgram(cl.loadSourceFromFile("ewald.cl"), replacements);
ewaldSumsKernel = cl::Kernel(program, "calculateEwaldCosSinSums"); ewaldSumsKernel = cl::Kernel(program, "calculateEwaldCosSinSums");
...@@ -809,7 +804,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -809,7 +804,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
int particle1, particle2; int particle1, particle2;
double chargeProd, sigma, epsilon; double chargeProd, sigma, epsilon;
force.getExceptionParameters(exceptions[i], particle1, particle2, chargeProd, sigma, epsilon); force.getExceptionParameters(exceptions[i], particle1, particle2, chargeProd, sigma, epsilon);
exceptionParamsVector[i] = (mm_float4) {(float) (138.935485*chargeProd), (float) sigma, (float) (4.0*epsilon), 0.0f}; exceptionParamsVector[i] = (mm_float4) {(float) (ONE_4PI_EPS0*chargeProd), (float) sigma, (float) (4.0*epsilon), 0.0f};
exceptionIndicesVector[i] = (mm_int4) {particle1, particle2, forceBufferCounter[particle1]++, forceBufferCounter[particle2]++}; exceptionIndicesVector[i] = (mm_int4) {particle1, particle2, forceBufferCounter[particle1]++, forceBufferCounter[particle2]++};
} }
exceptionParams->upload(exceptionParamsVector); exceptionParams->upload(exceptionParamsVector);
...@@ -1081,7 +1076,7 @@ void OpenCLCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOB ...@@ -1081,7 +1076,7 @@ void OpenCLCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOB
} }
posq.upload(); posq.upload();
params->upload(paramsVector); params->upload(paramsVector);
prefactor = -138.935485*((1.0/force.getSoluteDielectric())-(1.0/force.getSolventDielectric())); prefactor = -ONE_4PI_EPS0*((1.0/force.getSoluteDielectric())-(1.0/force.getSolventDielectric()));
bool useCutoff = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff); bool useCutoff = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff);
bool usePeriodic = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff && force.getNonbondedMethod() != GBSAOBCForce::CutoffNonPeriodic); bool usePeriodic = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff && force.getNonbondedMethod() != GBSAOBCForce::CutoffNonPeriodic);
string source = cl.loadSourceFromFile("gbsaObc2.cl"); string source = cl.loadSourceFromFile("gbsaObc2.cl");
......
platforms/opencl/src/kernels/coulombLennardJones.cl
View file @ 042ab1aa
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
if (r2 < CUTOFF_SQUARED) { if (r2 < CUTOFF_SQUARED) {
bool needCorrection = isExcluded && atom1 != atom2 && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS; bool needCorrection = isExcluded && atom1 != atom2 && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS;
if (!isExcluded || needCorrection) { if (!isExcluded || needCorrection) {
const float prefactor = 138.935485f*posq1.w*posq2.w*invR; const float prefactor = 138.935456f*posq1.w*posq2.w*invR;
float alphaR = EWALD_ALPHA*r; float alphaR = EWALD_ALPHA*r;
float erfcAlphaR = erfc(alphaR); float erfcAlphaR = erfc(alphaR);
float tempForce; float tempForce;
...@@ -47,11 +47,11 @@ if (!isExcluded) { ...@@ -47,11 +47,11 @@ if (!isExcluded) {
#endif #endif
#if HAS_COULOMB #if HAS_COULOMB
#ifdef USE_CUTOFF #ifdef USE_CUTOFF
const float prefactor = 138.935485f*posq1.w*posq2.w; const float prefactor = 138.935456f*posq1.w*posq2.w;
tempForce += prefactor*(invR - 2.0f*REACTION_FIELD_K*r2); tempForce += prefactor*(invR - 2.0f*REACTION_FIELD_K*r2);
tempEnergy += prefactor*(invR + REACTION_FIELD_K*r2 - REACTION_FIELD_C); tempEnergy += prefactor*(invR + REACTION_FIELD_K*r2 - REACTION_FIELD_C);
#else #else
const float prefactor = 138.935485f*posq1.w*posq2.w*invR; const float prefactor = 138.935456f*posq1.w*posq2.w*invR;
tempForce += prefactor; tempForce += prefactor;
tempEnergy += prefactor; tempEnergy += prefactor;
#endif #endif
......
platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp
View file @ 042ab1aa
...@@ -289,7 +289,7 @@ void testCoulombLennardJones() { ...@@ -289,7 +289,7 @@ void testCoulombLennardJones() {
customSystem.addParticle(1.0); customSystem.addParticle(1.0);
} }
NonbondedForce* standardNonbonded = new NonbondedForce(); NonbondedForce* standardNonbonded = new NonbondedForce();
CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935485*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)"); CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935456*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
customNonbonded->addPerParticleParameter("q"); customNonbonded->addPerParticleParameter("q");
customNonbonded->addPerParticleParameter("sigma"); customNonbonded->addPerParticleParameter("sigma");
customNonbonded->addPerParticleParameter("eps"); customNonbonded->addPerParticleParameter("eps");
......
platforms/opencl/tests/TestOpenCLNonbondedForce.cpp
View file @ 042ab1aa
...@@ -72,10 +72,10 @@ void testCoulomb() { ...@@ -72,10 +72,10 @@ void testCoulomb() {
context.setPositions(positions); context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces(); const vector<Vec3>& forces = state.getForces();
double force = 138.935485*(-0.75)/4.0; double force = ONE_4PI_EPS0*(-0.75)/4.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL); ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL); ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(138.935485*(-0.75)/2.0, state.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(-0.75)/2.0, state.getPotentialEnergy(), TOL);
} }
void testLJ() { void testLJ() {
...@@ -174,8 +174,8 @@ void testExclusionsAnd14() { ...@@ -174,8 +174,8 @@ void testExclusionsAnd14() {
context2.setPositions(positions); context2.setPositions(positions);
state = context2.getState(State::Forces | State::Energy); state = context2.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces(); const vector<Vec3>& forces2 = state.getForces();
force = 138.935485*4/(r*r); force = ONE_4PI_EPS0*4/(r*r);
energy = 138.935485*4/r; energy = ONE_4PI_EPS0*4/r;
if (i == 3) { if (i == 3) {
force /= 1.2; force /= 1.2;
energy /= 1.2; energy /= 1.2;
...@@ -217,13 +217,13 @@ void testCutoff() { ...@@ -217,13 +217,13 @@ void testCutoff() {
const vector<Vec3>& forces = state.getForces(); const vector<Vec3>& forces = state.getForces();
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0); const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0); const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force1 = 138.935485*(1.0)*(0.25-2.0*krf*2.0); const double force1 = ONE_4PI_EPS0*(1.0)*(0.25-2.0*krf*2.0);
const double force2 = 138.935485*(1.0)*(1.0-2.0*krf*1.0); const double force2 = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL); ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL); ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL); ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL);
const double energy1 = 138.935485*(1.0)*(0.5+krf*4.0-crf); const double energy1 = ONE_4PI_EPS0*(1.0)*(0.5+krf*4.0-crf);
const double energy2 = 138.935485*(1.0)*(1.0+krf*1.0-crf); const double energy2 = ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf);
ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL);
} }
...@@ -309,8 +309,8 @@ void testCutoff14() { ...@@ -309,8 +309,8 @@ void testCutoff14() {
const vector<Vec3>& forces2 = state.getForces(); const vector<Vec3>& forces2 = state.getForces();
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0); const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0); const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
force = 138.935485*q*q*(1.0/(r*r)-2.0*krf*r); force = ONE_4PI_EPS0*q*q*(1.0/(r*r)-2.0*krf*r);
energy = 138.935485*q*q*(1.0/r+krf*r*r-crf); energy = ONE_4PI_EPS0*q*q*(1.0/r+krf*r*r-crf);
if (i == 3) { if (i == 3) {
force /= 1.2; force /= 1.2;
energy /= 1.2; energy /= 1.2;
...@@ -353,11 +353,11 @@ void testPeriodic() { ...@@ -353,11 +353,11 @@ void testPeriodic() {
const double eps = 78.3; const double eps = 78.3;
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0); const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0); const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force = 138.935485*(1.0)*(1.0-2.0*krf*1.0); const double force = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL); ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL); ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL); ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
} }
......
platforms/reference/src/SimTKUtilities/SimTKOpenMMRealType.h
View file @ 042ab1aa
...@@ -138,8 +138,7 @@ ...@@ -138,8 +138,7 @@
#define ATOMICMASS_keV (940000.0) /* Atomic mass in keV */ #define ATOMICMASS_keV (940000.0) /* Atomic mass in keV */
#define ELECTRONMASS_keV (512.0) /* Electron mas in keV */ #define ELECTRONMASS_keV (512.0) /* Electron mas in keV */
#define FACEL 332.0636*CAL2JOULE /* (sqrt(ONE_4PI_EPS0)) */ #define ONE_4PI_EPS0 138.935456
#define ONE_4PI_EPS0 FACEL*0.1
#define PRESFAC (16.6054) /* bar / pressure unity */ #define PRESFAC (16.6054) /* bar / pressure unity */
#define ENM2DEBYE 48.0321 /* Convert electron nm to debye */ #define ENM2DEBYE 48.0321 /* Convert electron nm to debye */
#define DEBYE2ENM 0.02081941 #define DEBYE2ENM 0.02081941
......
platforms/reference/src/gbsa/CpuGBVI.cpp
View file @ 042ab1aa
...@@ -433,7 +433,6 @@ RealOpenMM CpuGBVI::computeBornEnergy( const RealOpenMM* bornRadii, RealOpenMM** ...@@ -433,7 +433,6 @@ RealOpenMM CpuGBVI::computeBornEnergy( const RealOpenMM* bornRadii, RealOpenMM**
static const RealOpenMM half = (RealOpenMM) 0.5; static const RealOpenMM half = (RealOpenMM) 0.5;
static const RealOpenMM fourth = (RealOpenMM) 0.25; static const RealOpenMM fourth = (RealOpenMM) 0.25;
static const RealOpenMM eighth = (RealOpenMM) 0.125; static const RealOpenMM eighth = (RealOpenMM) 0.125;
static const RealOpenMM CAL_TO_JOULE = (RealOpenMM) 0.4184;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -501,7 +500,7 @@ RealOpenMM e3 = -partialChargeI2*partialCharges[atomJ]*Sgb( t )/deltaR[Reference ...@@ -501,7 +500,7 @@ RealOpenMM e3 = -partialChargeI2*partialCharges[atomJ]*Sgb( t )/deltaR[Reference
energy += two*partialChargeI*atomIEnergy; energy += two*partialChargeI*atomIEnergy;
} }
energy *= 0.4184f*preFactor; energy *= preFactor;
energy -= cavityEnergy; energy -= cavityEnergy;
#if( GBVIDebug == 1 ) #if( GBVIDebug == 1 )
...@@ -551,7 +550,6 @@ int CpuGBVI::computeBornForces( const RealOpenMM* bornRadii, RealOpenMM** atomCo ...@@ -551,7 +550,6 @@ int CpuGBVI::computeBornForces( const RealOpenMM* bornRadii, RealOpenMM** atomCo
static const RealOpenMM oneThird = (RealOpenMM) (1.0/3.0); static const RealOpenMM oneThird = (RealOpenMM) (1.0/3.0);
static const RealOpenMM fourth = (RealOpenMM) 0.25; static const RealOpenMM fourth = (RealOpenMM) 0.25;
static const RealOpenMM eighth = (RealOpenMM) 0.125; static const RealOpenMM eighth = (RealOpenMM) 0.125;
static const RealOpenMM CAL_TO_JOULE = (RealOpenMM) 0.4184;
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -664,12 +662,12 @@ if( atomI == 0 ){ ...@@ -664,12 +662,12 @@ if( atomI == 0 ){
#if( GBVIDebug == 1 ) #if( GBVIDebug == 1 )
{ {
double stupidFactor = three/CAL_TO_JOULE; double stupidFactor = three;
RealOpenMM conversion = (RealOpenMM)(CAL_TO_JOULE*gbviParameters->getTau()); RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
int maxPrint = 10; int maxPrint = 10;
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii(); const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
(void) fprintf( logFile, "Conversion=%14.6e %14.6e*%14.6e (tau)\n", conversion, CAL_TO_JOULE, gbviParameters->getTau() ); (void) fprintf( logFile, "Conversion=%14.6e %14.6e*%14.6e (tau)\n", conversion, 1, gbviParameters->getTau() );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM R = atomicRadii[atomI]; RealOpenMM R = atomicRadii[atomI];
...@@ -717,7 +715,7 @@ if( atomI == 0 ){ ...@@ -717,7 +715,7 @@ if( atomI == 0 ){
#endif #endif
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii(); const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
RealOpenMM stupidFactor = three/0.4184f; RealOpenMM stupidFactor = three;
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM R = atomicRadii[atomI]; RealOpenMM R = atomicRadii[atomI];
...@@ -830,7 +828,7 @@ if( atomI == 0 ){ ...@@ -830,7 +828,7 @@ if( atomI == 0 ){
// convert from cal to Joule & apply prefactor tau = (1/diel_solute - 1/diel_solvent) // convert from cal to Joule & apply prefactor tau = (1/diel_solute - 1/diel_solvent)
RealOpenMM conversion = (RealOpenMM)(0.4184f*gbviParameters->getTau()); RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
inputForces[atomI][0] += conversion*forces[atomI][0]; inputForces[atomI][0] += conversion*forces[atomI][0];
inputForces[atomI][1] += conversion*forces[atomI][1]; inputForces[atomI][1] += conversion*forces[atomI][1];
......
platforms/reference/src/gbsa/CpuObc.cpp
View file @ 042ab1aa
...@@ -632,15 +632,12 @@ int CpuObc::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** atomCoo ...@@ -632,15 +632,12 @@ int CpuObc::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** atomCoo
// cal to Joule conversion // cal to Joule conversion
RealOpenMM conversion = (RealOpenMM)(0.4184);
// RealOpenMM conversion = (RealOpenMM)(0.4184/0.987);
//fprintf( stderr, "Obc includes .987 XXXXXXXXXXXXXXXXXXXXXXX\n" );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
inputForces[atomI][0] += conversion*forces[atomI][0]; inputForces[atomI][0] += forces[atomI][0];
inputForces[atomI][1] += conversion*forces[atomI][1]; inputForces[atomI][1] += forces[atomI][1];
inputForces[atomI][2] += conversion*forces[atomI][2]; inputForces[atomI][2] += forces[atomI][2];
} }
setEnergy( obcEnergy*conversion ); setEnergy( obcEnergy );
// copy new Born radii and obcChain values into permanent array // copy new Born radii and obcChain values into permanent array
...@@ -1259,7 +1256,7 @@ if( logFile && atomI >= 0 ){ ...@@ -1259,7 +1256,7 @@ if( logFile && atomI >= 0 ){
} }
RealOpenMM conversion = (RealOpenMM)0.4184; RealOpenMM conversion = (RealOpenMM) 1;
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
forces[atomI][0] *= conversion; forces[atomI][0] *= conversion;
forces[atomI][1] *= conversion; forces[atomI][1] *= conversion;
......
platforms/reference/src/gbsa/ImplicitSolventParameters.cpp
View file @ 042ab1aa
...@@ -637,12 +637,12 @@ void ImplicitSolventParameters::_initializeImplicitSolventConstants( void ){ ...@@ -637,12 +637,12 @@ void ImplicitSolventParameters::_initializeImplicitSolventConstants( void ){
_solventDielectric = (RealOpenMM) 78.3; _solventDielectric = (RealOpenMM) 78.3;
_kcalA_To_kJNm = (RealOpenMM) 0.4184; _kcalA_To_kJNm = (RealOpenMM) 0.4184;
_probeRadius = (RealOpenMM) 0.14; _probeRadius = (RealOpenMM) 0.14;
_electricConstant = (RealOpenMM) -166.02691; _electricConstant = (RealOpenMM) -0.5*ONE_4PI_EPS0;
//_pi4Asolv = (RealOpenMM) PI_M*4.0*0.0049*1000.0; //_pi4Asolv = (RealOpenMM) PI_M*4.0*0.0049*1000.0;
//_pi4Asolv = (RealOpenMM) PI_M*19.6; //_pi4Asolv = (RealOpenMM) PI_M*19.6;
//_pi4Asolv = (RealOpenMM) PI_M*4.0*0.0054; //_pi4Asolv = (RealOpenMM) PI_M*4.0*0.0054;
_pi4Asolv = (RealOpenMM) (PI_M*0.0216*1000.0); _pi4Asolv = (RealOpenMM) 28.3919551;
//_pi4Asolv = (RealOpenMM) -400.71504079; //_pi4Asolv = (RealOpenMM) -400.71504079;
//_pi4Asolv = (RealOpenMM) 0.0; //_pi4Asolv = (RealOpenMM) 0.0;
......
platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp
View file @ 042ab1aa
...@@ -250,7 +250,7 @@ void testCoulombLennardJones() { ...@@ -250,7 +250,7 @@ void testCoulombLennardJones() {
customSystem.addParticle(1.0); customSystem.addParticle(1.0);
} }
NonbondedForce* standardNonbonded = new NonbondedForce(); NonbondedForce* standardNonbonded = new NonbondedForce();
CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935485*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)"); CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935456*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
customNonbonded->addPerParticleParameter("q"); customNonbonded->addPerParticleParameter("q");
customNonbonded->addPerParticleParameter("sigma"); customNonbonded->addPerParticleParameter("sigma");
customNonbonded->addPerParticleParameter("eps"); customNonbonded->addPerParticleParameter("eps");
......
platforms/reference/tests/TestReferenceNonbondedForce.cpp
View file @ 042ab1aa
...@@ -66,10 +66,10 @@ void testCoulomb() { ...@@ -66,10 +66,10 @@ void testCoulomb() {
context.setPositions(positions); context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy); State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces(); const vector<Vec3>& forces = state.getForces();
double force = 138.935485*(-0.75)/4.0; double force = ONE_4PI_EPS0*(-0.75)/4.0;
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL); ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL); ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(138.935485*(-0.75)/2.0, state.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(-0.75)/2.0, state.getPotentialEnergy(), TOL);
} }
void testLJ() { void testLJ() {
...@@ -169,8 +169,8 @@ void testExclusionsAnd14() { ...@@ -169,8 +169,8 @@ void testExclusionsAnd14() {
context.setPositions(positions); context.setPositions(positions);
state = context.getState(State::Forces | State::Energy); state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces(); const vector<Vec3>& forces2 = state.getForces();
force = 138.935485*4/(r*r); force = ONE_4PI_EPS0*4/(r*r);
energy = 138.935485*4/r; energy = ONE_4PI_EPS0*4/r;
if (i == 3) { if (i == 3) {
force /= 1.2; force /= 1.2;
energy /= 1.2; energy /= 1.2;
...@@ -212,13 +212,13 @@ void testCutoff() { ...@@ -212,13 +212,13 @@ void testCutoff() {
const vector<Vec3>& forces = state.getForces(); const vector<Vec3>& forces = state.getForces();
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0); const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0); const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force1 = 138.935485*(1.0)*(0.25-2.0*krf*2.0); const double force1 = ONE_4PI_EPS0*(1.0)*(0.25-2.0*krf*2.0);
const double force2 = 138.935485*(1.0)*(1.0-2.0*krf*1.0); const double force2 = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL); ASSERT_EQUAL_VEC(Vec3(0, -force1, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL); ASSERT_EQUAL_VEC(Vec3(0, force1-force2, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL); ASSERT_EQUAL_VEC(Vec3(0, force2, 0), forces[2], TOL);
const double energy1 = 138.935485*(1.0)*(0.5+krf*4.0-crf); const double energy1 = ONE_4PI_EPS0*(1.0)*(0.5+krf*4.0-crf);
const double energy2 = 138.935485*(1.0)*(1.0+krf*1.0-crf); const double energy2 = ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf);
ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL);
} }
...@@ -304,8 +304,8 @@ void testCutoff14() { ...@@ -304,8 +304,8 @@ void testCutoff14() {
const vector<Vec3>& forces2 = state.getForces(); const vector<Vec3>& forces2 = state.getForces();
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0); const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0); const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
force = 138.935485*q*q*(1.0/(r*r)-2.0*krf*r); force = ONE_4PI_EPS0*q*q*(1.0/(r*r)-2.0*krf*r);
energy = 138.935485*q*q*(1.0/r+krf*r*r-crf); energy = ONE_4PI_EPS0*q*q*(1.0/r+krf*r*r-crf);
if (i == 3) { if (i == 3) {
force /= 1.2; force /= 1.2;
energy /= 1.2; energy /= 1.2;
...@@ -348,11 +348,11 @@ void testPeriodic() { ...@@ -348,11 +348,11 @@ void testPeriodic() {
const double eps = 78.3; const double eps = 78.3;
const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0); const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0); const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
const double force = 138.935485*(1.0)*(1.0-2.0*krf*1.0); const double force = ONE_4PI_EPS0*(1.0)*(1.0-2.0*krf*1.0);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL); ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL); ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL); ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
} }
int main() { int main() {
......
plugins/freeEnergy/platforms/reference/src/gbsa/CpuGBVISoftcore.cpp
View file @ 042ab1aa
...@@ -32,8 +32,6 @@ ...@@ -32,8 +32,6 @@
#include "../SimTKReference/ReferenceForce.h" #include "../SimTKReference/ReferenceForce.h"
#include <math.h> #include <math.h>
static const RealOpenMM CAL_TO_JOULE = 0.4184f;
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
CpuGBVISoftcore constructor CpuGBVISoftcore constructor
...@@ -688,7 +686,7 @@ RealOpenMM e3 = -partialChargeI2*partialCharges[atomJ]*Sgb( t )/deltaR[Reference ...@@ -688,7 +686,7 @@ RealOpenMM e3 = -partialChargeI2*partialCharges[atomJ]*Sgb( t )/deltaR[Reference
energy += two*partialChargeI*atomIEnergy; energy += two*partialChargeI*atomIEnergy;
} }
energy *= CAL_TO_JOULE*preFactor; energy *= preFactor;
energy -= cavityEnergy; energy -= cavityEnergy;
#if( GBVISoftcoreDebug == 1 ) #if( GBVISoftcoreDebug == 1 )
...@@ -850,13 +848,13 @@ if( atomI == 0 ){ ...@@ -850,13 +848,13 @@ if( atomI == 0 ){
#if( GBVISoftcoreDebug == 1 ) #if( GBVISoftcoreDebug == 1 )
{ {
double stupidFactor = three/CAL_TO_JOULE; double stupidFactor = three;
RealOpenMM conversion = (RealOpenMM)(CAL_TO_JOULE*gbviParameters->getTau()); RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
int maxPrint = 20; int maxPrint = 20;
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii(); const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
RealOpenMM* switchDeriviative = getSwitchDeriviative(); RealOpenMM* switchDeriviative = getSwitchDeriviative();
(void) fprintf( logFile, "F1: Conversion=%14.6e %14.6e*%14.6e (tau)\n", conversion, CAL_TO_JOULE, gbviParameters->getTau() ); (void) fprintf( logFile, "F1: Conversion=%14.6e %14.6e*%14.6e (tau)\n", conversion, 1, gbviParameters->getTau() );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM R = atomicRadii[atomI]; RealOpenMM R = atomicRadii[atomI];
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]); RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
...@@ -891,13 +889,13 @@ if( atomI == 0 ){ ...@@ -891,13 +889,13 @@ if( atomI == 0 ){
#if( GBVISoftcoreDebug == 2 ) #if( GBVISoftcoreDebug == 2 )
{ {
double stupidFactor = three/CAL_TO_JOULE; double stupidFactor = three;
RealOpenMM conversion = (RealOpenMM)(CAL_TO_JOULE*gbviParameters->getTau()); RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
int maxPrint = 1000000; int maxPrint = 1000000;
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii(); const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
RealOpenMM* switchDeriviative = getSwitchDeriviative(); RealOpenMM* switchDeriviative = getSwitchDeriviative();
(void) fprintf( logFile, "F1: Conversion=%14.6e %14.6e*%14.6e (tau)\n", conversion, CAL_TO_JOULE, gbviParameters->getTau() ); (void) fprintf( logFile, "F1: Conversion=%14.6e %14.6e*%14.6e (tau)\n", conversion, 1, gbviParameters->getTau() );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM R = atomicRadii[atomI]; RealOpenMM R = atomicRadii[atomI];
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]); RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
...@@ -943,7 +941,7 @@ if( atomI == 0 ){ ...@@ -943,7 +941,7 @@ if( atomI == 0 ){
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii(); const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
RealOpenMM* switchDeriviative = getSwitchDeriviative(); RealOpenMM* switchDeriviative = getSwitchDeriviative();
RealOpenMM stupidFactor = three/CAL_TO_JOULE; RealOpenMM stupidFactor = three;
const RealOpenMM* bornRadiusScaleFactors = gbviParameters->getBornRadiusScaleFactors(); const RealOpenMM* bornRadiusScaleFactors = gbviParameters->getBornRadiusScaleFactors();
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
...@@ -1054,7 +1052,7 @@ if( atomI == 0 ){ ...@@ -1054,7 +1052,7 @@ if( atomI == 0 ){
(void) fprintf( logFile, "\nPre conversion\n" ); (void) fprintf( logFile, "\nPre conversion\n" );
(void) fprintf( logFile, "Atom ScaledRadii BornRadii BornForce SwitchDrv Forces\n" ); (void) fprintf( logFile, "Atom ScaledRadii BornRadii BornForce SwitchDrv Forces\n" );
double forceSum[3] = { 0.0, 0.0, 0.0 }; double forceSum[3] = { 0.0, 0.0, 0.0 };
RealOpenMM conversion = (RealOpenMM)(CAL_TO_JOULE*gbviParameters->getTau()); RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii(); const RealOpenMM* scaledRadii = gbviParameters->getScaledRadii();
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
forceSum[0] += forces[atomI][0]; forceSum[0] += forces[atomI][0];
...@@ -1071,7 +1069,7 @@ if( atomI == 0 ){ ...@@ -1071,7 +1069,7 @@ if( atomI == 0 ){
// convert from cal to Joule & apply prefactor tau = (1/diel_solute - 1/diel_solvent) // convert from cal to Joule & apply prefactor tau = (1/diel_solute - 1/diel_solvent)
RealOpenMM conversion = (RealOpenMM)(CAL_TO_JOULE*gbviParameters->getTau()); RealOpenMM conversion = (RealOpenMM)(gbviParameters->getTau());
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
inputForces[atomI][0] += conversion*forces[atomI][0]; inputForces[atomI][0] += conversion*forces[atomI][0];
inputForces[atomI][1] += conversion*forces[atomI][1]; inputForces[atomI][1] += conversion*forces[atomI][1];
......
plugins/freeEnergy/platforms/reference/src/gbsa/CpuObcSoftcore.cpp
View file @ 042ab1aa
...@@ -699,13 +699,12 @@ int CpuObcSoftcore::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** ...@@ -699,13 +699,12 @@ int CpuObcSoftcore::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM**
// cal to Joule conversion // cal to Joule conversion
RealOpenMM conversion = (RealOpenMM)0.4184;
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){ for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
inputForces[atomI][0] += conversion*forces[atomI][0]; inputForces[atomI][0] += forces[atomI][0];
inputForces[atomI][1] += conversion*forces[atomI][1]; inputForces[atomI][1] += forces[atomI][1];
inputForces[atomI][2] += conversion*forces[atomI][2]; inputForces[atomI][2] += forces[atomI][2];
} }
setEnergy( obcEnergy*conversion ); setEnergy( obcEnergy );
#if 0 #if 0
{ {
...@@ -1347,13 +1346,7 @@ if( logFile && atomI >= 0 ){ ...@@ -1347,13 +1346,7 @@ if( logFile && atomI >= 0 ){
} }
RealOpenMM conversion = (RealOpenMM)0.4184; setEnergy( obcEnergy );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
forces[atomI][0] *= conversion;
forces[atomI][1] *= conversion;
forces[atomI][2] *= conversion;
}
setEnergy( obcEnergy*conversion );
if( 1 ){ if( 1 ){
......
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