Skip to content

GitLab

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

All platforms now handle nonbonded exceptions in a consistent way (they ignore... 

All platforms now handle nonbonded exceptions in a consistent way (they ignore cutoffs and periodic boundary conditions).  Also started OpenCL implementation of CustomHbondForce.
parent ca63e6d8
Show whitespace changes
Inline Side-by-side
Showing with 639 additions and 329 deletions
platforms/cuda/src/kernels/kCalculateCDLJForces.h
View file @ b97ff259
......@@ -213,7 +213,7 @@ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUnit)
dEdR *= invR * invR;
#ifdef USE_CUTOFF
#ifdef USE_EWALD
if ((!(excl & 0x1) && !needCorrection) || r2 > cSim.nonbondedCutoffSqr)
if (!needCorrection && (!(excl & 0x1) || r2 > cSim.nonbondedCutoffSqr))
#else
if (!(excl & 0x1) || r2 > cSim.nonbondedCutoffSqr)
#endif
......@@ -501,7 +501,7 @@ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUnit)
dEdR *= invR * invR;
#ifdef USE_CUTOFF
#ifdef USE_EWALD
if ((!(excl & 0x1) && !needCorrection) || r2 > cSim.nonbondedCutoffSqr)
if (!needCorrection && (!(excl & 0x1) || r2 > cSim.nonbondedCutoffSqr))
#else
if (!(excl & 0x1) || r2 > cSim.nonbondedCutoffSqr)
#endif
......
platforms/cuda/src/kernels/kCalculateLocalForces.cu
View file @ b97ff259
......@@ -461,8 +461,6 @@ void kCalculateLocalForces_kernel()
pos += blockDim.x * gridDim.x;
}
if (cSim.nonbondedMethod == NO_CUTOFF || cSim.nonbondedMethod == EWALD || cSim.nonbondedMethod == PARTICLE_MESH_EWALD)
{
while (pos < cSim.LJ14_offset)
{
unsigned int pos1 = pos - cSim.rb_dihedral_offset;
......@@ -506,121 +504,7 @@ void kCalculateLocalForces_kernel()
}
pos += blockDim.x * gridDim.x;
}
}
else if (cSim.nonbondedMethod == CUTOFF)
{
float LJ14_energy;
while (pos < cSim.LJ14_offset)
{
unsigned int pos1 = pos - cSim.rb_dihedral_offset;
if (pos1 < cSim.LJ14s)
{
int4 atom = cSim.pLJ14ID[pos1];
float4 LJ14 = cSim.pLJ14Parameter[pos1];
float4 a1 = cSim.pPosq[atom.x];
float4 a2 = cSim.pPosq[atom.y];
float3 d;
d.x = a1.x - a2.x;
d.y = a1.y - a2.y;
d.z = a1.z - a2.z;
float r2 = DOT3(d, d);
float inverseR = 1.0f / sqrt(r2);
float sig2 = inverseR * LJ14.y;
sig2 *= sig2;
float sig6 = sig2 * sig2 * sig2;
float dEdR = LJ14.x * (12.0f * sig6 - 6.0f) * sig6;
/* E */
LJ14_energy = LJ14.x * (sig6 - 1.0f) * sig6;
LJ14_energy += LJ14.z * (inverseR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
dEdR += LJ14.z * (inverseR - 2.0f * cSim.reactionFieldK * r2);
dEdR *= inverseR * inverseR;
if (r2 > cSim.nonbondedCutoffSqr)
{
dEdR = 0.0f;
/* E */
LJ14_energy = 0.0f;
}
/* E */
energy += LJ14_energy;
unsigned int offsetA = atom.x + atom.z * cSim.stride;
unsigned int offsetB = atom.y + atom.w * cSim.stride;
float4 forceA = cSim.pForce4[offsetA];
float4 forceB = cSim.pForce4[offsetB];
d.x *= dEdR;
d.y *= dEdR;
d.z *= dEdR;
forceA.x += d.x;
forceA.y += d.y;
forceA.z += d.z;
forceB.x -= d.x;
forceB.y -= d.y;
forceB.z -= d.z;
cSim.pForce4[offsetA] = forceA;
cSim.pForce4[offsetB] = forceB;
}
pos += blockDim.x * gridDim.x;
}
}
else if (cSim.nonbondedMethod == PERIODIC)
{
float LJ14_energy;
while (pos < cSim.LJ14_offset)
{
unsigned int pos1 = pos - cSim.rb_dihedral_offset;
if (pos1 < cSim.LJ14s)
{
int4 atom = cSim.pLJ14ID[pos1];
float4 LJ14 = cSim.pLJ14Parameter[pos1];
float4 a1 = cSim.pPosq[atom.x];
float4 a2 = cSim.pPosq[atom.y];
float3 d;
d.x = a1.x - a2.x;
d.y = a1.y - a2.y;
d.z = a1.z - a2.z;
d.x -= floor(d.x/cSim.periodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
d.y -= floor(d.y/cSim.periodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
d.z -= floor(d.z/cSim.periodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
float r2 = DOT3(d, d);
float inverseR = 1.0f / sqrt(r2);
float sig2 = inverseR * LJ14.y;
sig2 *= sig2;
float sig6 = sig2 * sig2 * sig2;
float dEdR = LJ14.x * (12.0f * sig6 - 6.0f) * sig6;
/* E */
LJ14_energy = LJ14.x * (sig6 - 1.0f) * sig6;
LJ14_energy += LJ14.z * (inverseR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
dEdR += LJ14.z * (inverseR - 2.0f * cSim.reactionFieldK * r2);
dEdR *= inverseR * inverseR;
if (r2 > cSim.nonbondedCutoffSqr)
{
dEdR = 0.0f;
/* E */
LJ14_energy = 0.0f;
}
/* E */
energy += LJ14_energy;
unsigned int offsetA = atom.x + atom.z * cSim.stride;
unsigned int offsetB = atom.y + atom.w * cSim.stride;
float4 forceA = cSim.pForce4[offsetA];
float4 forceB = cSim.pForce4[offsetB];
d.x *= dEdR;
d.y *= dEdR;
d.z *= dEdR;
forceA.x += d.x;
forceA.y += d.y;
forceA.z += d.z;
forceB.x -= d.x;
forceB.y -= d.y;
forceB.z -= d.z;
cSim.pForce4[offsetA] = forceA;
cSim.pForce4[offsetB] = forceB;
}
pos += blockDim.x * gridDim.x;
}
}
cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += energy;
}
......
platforms/cuda/tests/TestCudaEwald.cpp
View file @ b97ff259
......@@ -79,17 +79,20 @@ void testEwaldPME(bool includeExceptions) {
nonbonded->addParticle(1.0, 1.0,0.0);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(-1.0, 1.0,0.0);
if (includeExceptions) {
// Add some exclusions.
for (int i = 0; i < numParticles-1; i++)
nonbonded->addException(i, i+1, 0.0, 1.0, 0.0);
}
system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
vector<Vec3> positions(numParticles);
#include "nacl_amorph.dat"
if (includeExceptions) {
// Add some exclusions.
for (int i = 0; i < numParticles-1; i++) {
Vec3 delta = positions[i]-positions[i+1];
if (sqrt(delta.dot(delta)) < 0.5*cutoff)
nonbonded->addException(i, i+1, i%2 == 0 ? 0.0 : 0.5, 1.0, 0.0);
}
}
// (1) Check whether the Reference and Cuda platforms agree when using Ewald Method
......
platforms/cuda/tests/TestCudaNonbondedForce.cpp
View file @ b97ff259
......@@ -312,10 +312,8 @@ void testCutoff14() {
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces();
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);
force = ONE_4PI_EPS0*q*q*(1.0/(r*r)-2.0*krf*r);
energy = ONE_4PI_EPS0*q*q*(1.0/r+krf*r*r-crf);
force = ONE_4PI_EPS0*q*q/(r*r);
energy = ONE_4PI_EPS0*q*q/r;
if (i == 3) {
force /= 1.2;
energy /= 1.2;
......
platforms/opencl/src/OpenCLKernelFactory.cpp
View file @ b97ff259
......@@ -62,6 +62,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
return new OpenCLCalcCustomGBForceKernel(name, platform, cl, context.getSystem());
if (name == CalcCustomExternalForceKernel::Name())
return new OpenCLCalcCustomExternalForceKernel(name, platform, cl, context.getSystem());
if (name == CalcCustomHbondForceKernel::Name())
return new OpenCLCalcCustomHbondForceKernel(name, platform, cl, context.getSystem());
if (name == IntegrateVerletStepKernel::Name())
return new OpenCLIntegrateVerletStepKernel(name, platform, cl);
if (name == IntegrateLangevinStepKernel::Name())
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ b97ff259
......@@ -37,6 +37,7 @@
#include "lepton/Parser.h"
#include "lepton/ParsedExpression.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "lepton/Operation.h"
#include <cmath>
using namespace OpenMM;
......@@ -55,6 +56,13 @@ static string intToString(int value) {
return s.str();
}
static bool isZeroExpression(const Lepton::ParsedExpression& expression) {
const Lepton::Operation& op = expression.getRootNode().getOperation();
if (op.getId() != Lepton::Operation::CONSTANT)
return false;
return (dynamic_cast<const Lepton::Operation::Constant&>(op).getValue() == 0.0);
}
void OpenCLCalcForcesAndEnergyKernel::initialize(const System& system) {
}
......@@ -318,9 +326,6 @@ void OpenCLCalcCustomBondForceKernel::initialize(const System& system, const Cus
// Record information for the expressions.
vector<string> paramNames;
for (int i = 0; i < force.getNumPerBondParameters(); i++)
paramNames.push_back(force.getPerBondParameterName(i));
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
......@@ -556,9 +561,6 @@ void OpenCLCalcCustomAngleForceKernel::initialize(const System& system, const Cu
// Record information for the expressions.
vector<string> paramNames;
for (int i = 0; i < force.getNumPerAngleParameters(); i++)
paramNames.push_back(force.getPerAngleParameterName(i));
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
......@@ -880,9 +882,6 @@ void OpenCLCalcCustomTorsionForceKernel::initialize(const System& system, const
// Record information for the expressions.
vector<string> paramNames;
for (int i = 0; i < force.getNumPerTorsionParameters(); i++)
paramNames.push_back(force.getPerTorsionParameterName(i));
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
......@@ -1218,7 +1217,6 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source);
if (hasLJ)
cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo("sigmaEpsilon", "float2", sizeof(cl_float2), sigmaEpsilon->getDeviceBuffer()));
cutoffSquared = force.getCutoffDistance()*force.getCutoffDistance();
// Initialize the exceptions.
......@@ -1243,11 +1241,9 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
maxBuffers = max(maxBuffers, forceBufferCounter[i]);
}
cl.addForce(new OpenCLNonbondedForceInfo(maxBuffers, force));
if (useCutoff) {
defines["USE_CUTOFF"] = "1";
}
if (usePeriodic)
defines["USE_PERIODIC"] = "1";
defines.clear();
defines["NUM_ATOMS"] = intToString(numParticles);
defines["NUM_EXCEPTIONS"] = intToString(numExceptions);
cl::Program program = cl.createProgram(OpenCLKernelSources::nonbondedExceptions, defines);
exceptionsKernel = cl::Kernel(program, "computeNonbondedExceptions");
}
......@@ -1256,16 +1252,11 @@ void OpenCLCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
if (!hasInitializedKernel) {
hasInitializedKernel = true;
if (exceptionIndices != NULL) {
int numExceptions = exceptionIndices->getSize();
exceptionsKernel.setArg<cl_int>(0, cl.getPaddedNumAtoms());
exceptionsKernel.setArg<cl_int>(1, numExceptions);
exceptionsKernel.setArg<cl_float>(2, (cl_float) cutoffSquared);
exceptionsKernel.setArg<mm_float4>(3, cl.getNonbondedUtilities().getPeriodicBoxSize());
exceptionsKernel.setArg<cl::Buffer>(4, cl.getForceBuffers().getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(5, cl.getEnergyBuffer().getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(6, cl.getPosq().getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(7, exceptionParams->getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(8, exceptionIndices->getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(0, cl.getForceBuffers().getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(1, cl.getEnergyBuffer().getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(2, cl.getPosq().getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(3, exceptionParams->getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(4, exceptionIndices->getDeviceBuffer());
}
if (cosSinSums != NULL) {
ewaldSumsKernel.setArg<cl::Buffer>(0, cl.getEnergyBuffer().getDeviceBuffer());
......@@ -1438,9 +1429,6 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
// Record information for the expressions.
vector<string> paramNames;
for (int i = 0; i < force.getNumPerParticleParameters(); i++)
paramNames.push_back(force.getPerParticleParameterName(i));
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
......@@ -1782,9 +1770,6 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
// Record the global parameters.
vector<string> paramNames;
for (int i = 0; i < force.getNumPerParticleParameters(); i++)
paramNames.push_back(force.getPerParticleParameterName(i));
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
......@@ -2375,9 +2360,6 @@ void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const
// Record information for the expressions.
vector<string> paramNames;
for (int i = 0; i < force.getNumPerParticleParameters(); i++)
paramNames.push_back(force.getPerParticleParameterName(i));
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
......@@ -2461,6 +2443,306 @@ double OpenCLCalcCustomExternalForceKernel::executeEnergy(ContextImpl& context)
return 0.0;
}
class OpenCLCustomHbondForceInfo : public OpenCLForceInfo {
public:
OpenCLCustomHbondForceInfo(int requiredBuffers, const CustomHbondForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
}
bool areParticlesIdentical(int particle1, int particle2) {
return true;
}
int getNumParticleGroups() {
return force.getNumDonors()+force.getNumAcceptors()+force.getNumExclusions();
}
void getParticlesInGroup(int index, std::vector<int>& particles) {
int p1, p2, p3;
vector<double> parameters;
if (index < force.getNumDonors()) {
force.getDonorParameters(index, p1, p2, p3, parameters);
particles.resize(3);
particles[0] = p1;
particles[1] = p2;
particles[2] = p3;
return;
}
index -= force.getNumDonors();
if (index < force.getNumAcceptors()) {
force.getAcceptorParameters(index, p1, p2, p3, parameters);
particles.resize(3);
particles[0] = p1;
particles[1] = p2;
particles[2] = p3;
return;
}
index -= force.getNumAcceptors();
int donor, acceptor;
force.getExclusionParticles(index, donor, acceptor);
particles.resize(6);
force.getDonorParameters(donor, p1, p2, p3, parameters);
particles[0] = p1;
particles[1] = p2;
particles[2] = p3;
force.getAcceptorParameters(acceptor, p1, p2, p3, parameters);
particles[3] = p1;
particles[4] = p2;
particles[5] = p3;
}
bool areGroupsIdentical(int group1, int group2) {
int p1, p2, p3;
vector<double> params1, params2;
if (group1 < force.getNumDonors() && group2 < force.getNumDonors()) {
force.getDonorParameters(group1, p1, p2, p3, params1);
force.getDonorParameters(group2, p1, p2, p3, params2);
return (params1 == params2 && params1 == params2);
}
if (group1 < force.getNumDonors() || group2 < force.getNumDonors())
return false;
group1 -= force.getNumDonors();
group2 -= force.getNumDonors();
if (group1 < force.getNumAcceptors() && group2 < force.getNumAcceptors()) {
force.getAcceptorParameters(group1, p1, p2, p3, params1);
force.getAcceptorParameters(group2, p1, p2, p3, params2);
return (params1 == params2 && params1 == params2);
}
if (group1 < force.getNumAcceptors() || group2 < force.getNumAcceptors())
return false;
return true;
}
private:
const CustomHbondForce& force;
};
OpenCLCalcCustomHbondForceKernel::~OpenCLCalcCustomHbondForceKernel() {
if (donorParams != NULL)
delete donorParams;
if (acceptorParams != NULL)
delete acceptorParams;
if (donors != NULL)
delete donors;
if (acceptors != NULL)
delete acceptors;
if (donorBufferIndices != NULL)
delete donorBufferIndices;
if (acceptorBufferIndices != NULL)
delete acceptorBufferIndices;
if (globals != NULL)
delete globals;
if (tabulatedFunctionParams != NULL)
delete tabulatedFunctionParams;
for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
delete tabulatedFunctions[i];
}
void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const CustomHbondForce& force) {
int forceIndex;
for (forceIndex = 0; forceIndex < system.getNumForces() && &system.getForce(forceIndex) != &force; ++forceIndex)
;
string prefix = "custom"+intToString(forceIndex)+"_";
// Record the lists of donors and acceptors, and the parameters for each one.
numDonors = force.getNumDonors();
numAcceptors = force.getNumAcceptors();
int numParticles = system.getNumParticles();
donors = new OpenCLArray<mm_int4>(cl, numDonors, "customHbondDonors");
acceptors = new OpenCLArray<mm_int4>(cl, numAcceptors, "customHbondAcceptors");
donorParams = new OpenCLParameterSet(cl, force.getNumPerDonorParameters(), numDonors, "customHbondDonorParameters");
acceptorParams = new OpenCLParameterSet(cl, force.getNumPerAcceptorParameters(), numAcceptors, "customHbondAcceptorParameters");
if (force.getNumGlobalParameters() > 0)
globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customHbondGlobals", false, CL_MEM_READ_ONLY);
vector<vector<cl_float> > donorParamVector(numDonors);
vector<mm_int4> donorVector(numDonors);
for (int i = 0; i < numDonors; i++) {
vector<double> parameters;
force.getDonorParameters(i, donorVector[i].x, donorVector[i].y, donorVector[i].z, parameters);
donorParamVector[i].resize(parameters.size());
for (int j = 0; j < (int) parameters.size(); j++)
donorParamVector[i][j] = (cl_float) parameters[j];
}
donors->upload(donorVector);
donorParams->setParameterValues(donorParamVector);
vector<vector<cl_float> > acceptorParamVector(numAcceptors);
vector<mm_int4> acceptorVector(numAcceptors);
for (int i = 0; i < numAcceptors; i++) {
vector<double> parameters;
force.getAcceptorParameters(i, acceptorVector[i].x, acceptorVector[i].y, acceptorVector[i].z, parameters);
acceptorParamVector[i].resize(parameters.size());
for (int j = 0; j < (int) parameters.size(); j++)
acceptorParamVector[i][j] = (cl_float) parameters[j];
}
acceptors->upload(acceptorVector);
acceptorParams->setParameterValues(acceptorParamVector);
// Select a output buffer indices for each donor and acceptor.
donorBufferIndices = new OpenCLArray<mm_int4>(cl, numDonors, "customHbondDonorBuffers");
acceptorBufferIndices = new OpenCLArray<mm_int4>(cl, numAcceptors, "customHbondAcceptorBuffers");
vector<mm_int4> donorBufferVector(numDonors);
vector<mm_int4> acceptorBufferVector(numAcceptors);
vector<int> particleBuffers(numParticles, 0);
for (int i = 0; i < numDonors; i++)
donorBufferVector[i] = mm_int4(particleBuffers[donorVector[i].x]++, particleBuffers[donorVector[i].y]++, particleBuffers[donorVector[i].z]++, 0);
for (int i = 0; i < numAcceptors; i++)
acceptorBufferVector[i] = mm_int4(particleBuffers[acceptorVector[i].x]++, particleBuffers[acceptorVector[i].y]++, particleBuffers[acceptorVector[i].z]++, 0);
donorBufferIndices->upload(donorBufferVector);
acceptorBufferIndices->upload(acceptorBufferVector);
// Record exclusions.
vector<vector<int> > exclusionList(numDonors);
for (int i = 0; i < force.getNumExclusions(); i++) {
int donor, acceptor;
force.getExclusionParticles(i, donor, acceptor);
exclusionList[donor].push_back(acceptor);
}
// Record the tabulated functions.
OpenCLExpressionUtilities::FunctionPlaceholder fp;
map<string, Lepton::CustomFunction*> functions;
vector<pair<string, string> > functionDefinitions;
vector<mm_float4> tabulatedFunctionParamsVec(force.getNumFunctions());
stringstream tableArgs;
for (int i = 0; i < force.getNumFunctions(); i++) {
string name;
vector<double> values;
double min, max;
bool interpolating;
force.getFunctionParameters(i, name, values, min, max, interpolating);
string arrayName = prefix+"table"+intToString(i);
functionDefinitions.push_back(make_pair(name, arrayName));
functions[name] = &fp;
tabulatedFunctionParamsVec[i] = mm_float4((float) min, (float) max, (float) ((values.size()-1)/(max-min)), 0.0f);
vector<mm_float4> f = OpenCLExpressionUtilities::computeFunctionCoefficients(values, interpolating);
tabulatedFunctions.push_back(new OpenCLArray<mm_float4>(cl, values.size()-1, "TabulatedFunction"));
tabulatedFunctions[tabulatedFunctions.size()-1]->upload(f);
tableArgs << ", __global float4* " << arrayName;
}
if (force.getNumFunctions() > 0) {
tabulatedFunctionParams = new OpenCLArray<mm_float4>(cl, tabulatedFunctionParamsVec.size(), "tabulatedFunctionParameters", false, CL_MEM_READ_ONLY);
tabulatedFunctionParams->upload(tabulatedFunctionParamsVec);
tableArgs << ", __constant float4* " << prefix << "functionParams";
}
// Record information for the expressions.
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
globalParamNames[i] = force.getGlobalParameterName(i);
globalParamValues[i] = (cl_float) force.getGlobalParameterDefaultValue(i);
}
if (globals != NULL)
globals->upload(globalParamValues);
bool useCutoff = (force.getNonbondedMethod() != CustomHbondForce::NoCutoff);
bool usePeriodic = (force.getNonbondedMethod() != CustomHbondForce::NoCutoff && force.getNonbondedMethod() != CustomHbondForce::CutoffNonPeriodic);
Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
Lepton::ParsedExpression rForceExpression = energyExpression.differentiate("r").optimize();
Lepton::ParsedExpression thetaForceExpression = energyExpression.differentiate("theta").optimize();
Lepton::ParsedExpression psiForceExpression = energyExpression.differentiate("psi").optimize();
Lepton::ParsedExpression chiForceExpression = energyExpression.differentiate("chi").optimize();
map<string, Lepton::ParsedExpression> forceExpressions;
forceExpressions["energy += "] = energyExpression;
forceExpressions["float dEdR = "] = rForceExpression;
forceExpressions["float dEdTheta = "] = thetaForceExpression;
forceExpressions["float dEdPsi = "] = psiForceExpression;
forceExpressions["float dEdChi = "] = chiForceExpression;
// Create the kernels.
map<string, string> variables;
variables["r"] = "r";
variables["theta"] = "theta";
variables["psi"] = "psi";
variables["chi"] = "chi";
for (int i = 0; i < force.getNumPerDonorParameters(); i++) {
const string& name = force.getPerDonorParameterName(i);
variables[name] = "donorParams"+donorParams->getParameterSuffix(i);
}
for (int i = 0; i < force.getNumPerAcceptorParameters(); i++) {
const string& name = force.getPerAcceptorParameterName(i);
variables[name] = "acceptorParams"+acceptorParams->getParameterSuffix(i);
}
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
variables[name] = "globals["+intToString(i)+"]";
}
stringstream compute, extraArgs;
if (force.getNumGlobalParameters() > 0)
extraArgs << ", __constant float* globals";
for (int i = 0; i < (int) donorParams->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = donorParams->getBuffers()[i];
extraArgs << ", __global "+buffer.getType()+"* donor"+buffer.getName();
compute<<buffer.getType()<<" donorParams"<<(i+1)<<" = donor"<<buffer.getName()<<"[index];\n";
}
for (int i = 0; i < (int) acceptorParams->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = acceptorParams->getBuffers()[i];
extraArgs << ", __global "+buffer.getType()+"* acceptor"+buffer.getName();
compute<<buffer.getType()<<" acceptorParams"<<(i+1)<<" = acceptor"<<buffer.getName()<<"[index];\n";
}
compute << OpenCLExpressionUtilities::createExpressions(forceExpressions, variables, functionDefinitions, "temp", "functionParams");
map<string, string> replacements;
replacements["COMPUTE_FORCE"] = compute.str();
replacements["PARAMETER_ARGUMENTS"] = extraArgs.str()+tableArgs.str();
map<string, string> defines;
defines["PADDED_NUM_ATOMS"] = intToString(cl.getPaddedNumAtoms());
defines["NUM_DONORS"] = intToString(force.getNumDonors());
defines["NUM_ACCEPTORS"] = intToString(force.getNumAcceptors());
defines["M_PI"] = doubleToString(M_PI);
if (!isZeroExpression(rForceExpression))
defines["INCLUDE_R"] = "1";
if (!isZeroExpression(thetaForceExpression))
defines["INCLUDE_THETA"] = "1";
if (!isZeroExpression(psiForceExpression))
defines["INCLUDE_PSI"] = "1";
if (!isZeroExpression(chiForceExpression))
defines["INCLUDE_CHI"] = "1";
cl::Program program = cl.createProgram(cl.replaceStrings(OpenCLKernelSources::customHbondForce, replacements), defines);
kernel = cl::Kernel(program, "computeHbonds");
}
void OpenCLCalcCustomHbondForceKernel::executeForces(ContextImpl& context) {
if (globals != NULL) {
bool changed = false;
for (int i = 0; i < (int) globalParamNames.size(); i++) {
cl_float value = (cl_float) context.getParameter(globalParamNames[i]);
if (value != globalParamValues[i])
changed = true;
globalParamValues[i] = value;
}
if (changed)
globals->upload(globalParamValues);
}
if (!hasInitializedKernel) {
hasInitializedKernel = true;
int index = 0;
kernel.setArg<cl::Buffer>(index++, cl.getForceBuffers().getDeviceBuffer());
kernel.setArg<cl::Buffer>(index++, cl.getEnergyBuffer().getDeviceBuffer());
kernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
kernel.setArg<cl::Buffer>(index++, donors->getDeviceBuffer());
kernel.setArg<cl::Buffer>(index++, acceptors->getDeviceBuffer());
kernel.setArg<cl::Buffer>(index++, donorBufferIndices->getDeviceBuffer());
kernel.setArg<cl::Buffer>(index++, acceptorBufferIndices->getDeviceBuffer());
kernel.setArg(index++, OpenCLContext::ThreadBlockSize*sizeof(mm_float4), NULL);
kernel.setArg(index++, OpenCLContext::ThreadBlockSize*sizeof(mm_float4), NULL);
if (globals != NULL)
kernel.setArg<cl::Buffer>(index++, globals->getDeviceBuffer());
for (int i = 0; i < (int) donorParams->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = donorParams->getBuffers()[i];
kernel.setArg<cl::Buffer>(index++, buffer.getBuffer());
}
for (int i = 0; i < (int) acceptorParams->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = acceptorParams->getBuffers()[i];
kernel.setArg<cl::Buffer>(index++, buffer.getBuffer());
}
}
cl.executeKernel(kernel, std::max(numDonors, numAcceptors));
}
double OpenCLCalcCustomHbondForceKernel::executeEnergy(ContextImpl& context) {
executeForces(context);
return 0.0;
}
OpenCLIntegrateVerletStepKernel::~OpenCLIntegrateVerletStepKernel() {
}
......
platforms/opencl/src/OpenCLKernels.h
View file @ b97ff259
......@@ -488,7 +488,7 @@ private:
cl::Kernel pmeConvolutionKernel;
cl::Kernel pmeInterpolateForceKernel;
std::map<std::string, std::string> pmeDefines;
double cutoffSquared, ewaldSelfEnergy;
double ewaldSelfEnergy;
static const int PmeOrder = 5;
};
......@@ -664,6 +664,55 @@ private:
cl::Kernel kernel;
};
/**
* This kernel is invoked by CustomHbondForce to calculate the forces acting on the system.
*/
class OpenCLCalcCustomHbondForceKernel : public CalcCustomHbondForceKernel {
public:
OpenCLCalcCustomHbondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomHbondForceKernel(name, platform),
hasInitializedKernel(false), cl(cl), donorParams(NULL), acceptorParams(NULL), donors(NULL), acceptors(NULL),
donorBufferIndices(NULL), acceptorBufferIndices(NULL), globals(NULL), tabulatedFunctionParams(NULL), system(system) {
}
~OpenCLCalcCustomHbondForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the CustomHbondForce this kernel will be used for
*/
void initialize(const System& system, const CustomHbondForce& force);
/**
* Execute the kernel to calculate the forces.
*
* @param context the context in which to execute this kernel
*/
void executeForces(ContextImpl& context);
/**
* Execute the kernel to calculate the energy.
*
* @param context the context in which to execute this kernel
* @return the potential energy due to the CustomHbondForce
*/
double executeEnergy(ContextImpl& context);
private:
int numDonors, numAcceptors;
bool hasInitializedKernel;
OpenCLContext& cl;
OpenCLParameterSet* donorParams;
OpenCLParameterSet* acceptorParams;
OpenCLArray<cl_float>* globals;
OpenCLArray<mm_int4>* donors;
OpenCLArray<mm_int4>* acceptors;
OpenCLArray<mm_int4>* donorBufferIndices;
OpenCLArray<mm_int4>* acceptorBufferIndices;
OpenCLArray<mm_float4>* tabulatedFunctionParams;
std::vector<std::string> globalParamNames;
std::vector<cl_float> globalParamValues;
std::vector<OpenCLArray<mm_float4>*> tabulatedFunctions;
System& system;
cl::Kernel kernel;
};
/**
* This kernel is invoked by VerletIntegrator to take one time step.
*/
......
platforms/opencl/src/OpenCLPlatform.cpp
View file @ b97ff259
......@@ -58,6 +58,7 @@ OpenCLPlatform::OpenCLPlatform() {
registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
registerKernelFactory(CalcCustomGBForceKernel::Name(), factory);
registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
......
platforms/opencl/src/kernels/coulombLennardJones.cl
View file @ b97ff259
#if USE_EWALD
if (r2 < CUTOFF_SQUARED) {
bool needCorrection = isExcluded && atom1 != atom2 && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS;
if (!isExcluded || needCorrection) {
bool needCorrection = isExcluded && atom1 != atom2 && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS;
if (!isExcluded || needCorrection) {
const float prefactor = 138.935456f*posq1.w*posq2.w*invR;
float alphaR = EWALD_ALPHA*r;
float erfcAlphaR = erfc(alphaR);
float tempForce;
float tempForce = 0.0f;
if (needCorrection) {
// Subtract off the part of this interaction that was included in the reciprocal space contribution.
tempForce = -prefactor*((1.0f-erfcAlphaR)-alphaR*exp(-alphaR*alphaR)*TWO_OVER_SQRT_PI);
tempEnergy += -prefactor*(1.0f-erfcAlphaR);
}
else {
else if (r2 < CUTOFF_SQUARED) {
#if HAS_LENNARD_JONES
float sig = sigmaEpsilon1.x + sigmaEpsilon2.x;
float sig2 = invR*sig;
......@@ -27,7 +26,6 @@ if (r2 < CUTOFF_SQUARED) {
#endif
}
dEdR += tempForce*invR*invR;
}
}
#else
#ifdef USE_CUTOFF
......
platforms/opencl/src/kernels/customHbondForce.cl 0 → 100644
View file @ b97ff259
/**
* Compute the difference between two vectors, setting the fourth component to the squared magnitude.
*/
float4 delta(float4 vec1, float4 vec2) {
float4 result = (float4) (vec1.x-vec2.x, vec1.y-vec2.y, vec1.z-vec2.z, 0.0f);
result.w = result.x*result.x + result.y*result.y + result.z*result.z;
return result;
}
/**
* Compute the difference between two vectors, taking periodic boundary conditions into account
* and setting the fourth component to the squared magnitude.
*/
float4 deltaPeriodic(float4 vec1, float4 vec2) {
float4 result = (float4) (vec1.x-vec2.x, vec1.y-vec2.y, vec1.z-vec2.z, 0.0f);
#ifdef USE_PERIODIC
result.x -= floor(result.x/PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
result.y -= floor(result.y/PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
result.z -= floor(result.z/PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
#endif
result.w = result.x*result.x + result.y*result.y + result.z*result.z;
return result;
}
/**
* Compute the angle between two vectors. The w component of each vector should contain the squared magnitude.
*/
float computeAngle(float4 vec1, float4 vec2) {
float dot = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
float cosine = dot/sqrt(vec1.w*vec2.w);
float angle;
if (cosine > 0.99f || cosine < -0.99f) {
// We're close to the singularity in acos(), so take the cross product and use asin() instead.
float4 cross_prod = cross(vec1, vec2);
float scale = vec1.w*vec2.w;
angle = asin(sqrt(dot(cross_prod, cross_prod)/scale));
if (cosine < 0.0f)
angle = M_PI-angle;
}
else
angle = acos(cosine);
return angle;
}
/**
* Compute hbond interactions.
*/
__kernel void computeHbonds(__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, /*__global unsigned int* exclusions,
__global unsigned int* exclusionIndices, */__global int4* donorAtoms, __global int4* acceptorAtoms, __global int4* donorBufferIndices, __global int4* acceptorBufferIndices, __local float4* posBuffer, __local float4* deltaBuffer
PARAMETER_ARGUMENTS) {
float energy = 0.0f;
unsigned int tgx = get_local_id(0) & (get_local_size(0)-1);
unsigned int tbx = get_local_id(0) - tgx;
float4 f1 = 0;
float4 f2 = 0;
for (int donorIndex = get_global_id(0); donorIndex < NUM_DONORS; donorIndex += get_global_size(0)) {
// Load information about the donor this thread will compute forces on.
int4 atoms = donorAtoms[donorIndex];
float4 d1 = posq[atoms.x];
float4 d2 = posq[atoms.y];
float4 d3 = posq[atoms.z];
float4 deltaD1D2 = delta(d1, d2);
for (int acceptorStart = 0; acceptorStart < NUM_ACCEPTORS; acceptorStart += get_local_size(0)) {
// Load the next block of acceptors into local memory.
int blockSize = min((int) get_local_size(0), NUM_ACCEPTORS-acceptorStart);
if (tgx < blockSize) {
int4 atoms2 = acceptorAtoms[acceptorStart+tgx];
float4 pos1 = posq[atoms2.x];
float4 pos2 = posq[atoms2.y];
float4 pos3 = posq[atoms2.z];
posBuffer[get_local_id(0)] = pos1;
deltaBuffer[get_local_id(0)] = delta(pos2, pos1);
}
barrier(CLK_LOCAL_MEM_FENCE);
for (int index = 0; index < blockSize; index++) {
// Compute the interaction between a donor and an acceptor.
float4 a1 = posBuffer[index];
float4 deltaD1A1 = deltaPeriodic(d1, a1);
#ifdef USE_CUTOFF
if (deltaD1A1.w < CUTOFF_SQUARED) {
#endif
// Compute variables the force can depend on.
float r = sqrt(deltaD1A1.w);
float4 deltaA2A1 = deltaBuffer[index];
float theta = computeAngle(deltaD1A1, deltaD1D2);
float psi = computeAngle(deltaD1A1, deltaA2A1);
float4 cross1 = cross(deltaA2A1, deltaD1A1);
float4 cross2 = cross(deltaD1A1, deltaD1D2);
cross1.w = cross1.x*cross1.x + cross1.y*cross1.y + cross1.z*cross1.z;
cross2.w = cross2.x*cross2.x + cross2.y*cross2.y + cross2.z*cross2.z;
float chi = computeAngle(cross1, cross2);
chi = (dot(deltaA2A1, cross2) < 0 ? -chi : chi);
COMPUTE_FORCE
#ifdef INCLUDE_R
// Apply forces based on r.
f1.xyz -= (dEdR/r)*deltaD1A1.xyz;
#endif
#ifdef INCLUDE_THETA
// Apply forces based on theta.
float4 thetaCross = cross(deltaD1D2, deltaD1A1);
float lengthThetaCross = max(length(thetaCross), 1e-6f);
float4 deltaCross0 = cross(deltaD1D2, thetaCross)*dEdTheta/(deltaD1D2.w*lengthThetaCross);
float4 deltaCross2 = -cross(deltaD1A1, thetaCross)*dEdTheta/(deltaD1A1.w*lengthThetaCross);
float4 deltaCross1 = -(deltaCross0+deltaCross2);
f1.xyz += deltaCross1.xyz;
f2.xyz += deltaCross0.xyz;
#endif
#ifdef INCLUDE_PSI
// Apply forces based on psi.
float4 psiCross = cross(deltaA2A1, deltaD1A1);
float lengthPsiCross = max(length(psiCross), 1e-6f);
deltaCross0 = cross(deltaD1A1, psiCross)*dEdPsi/(deltaD1A1.w*lengthPsiCross);
// float4 deltaCross2 = -cross(deltaA2A1, psiCross)*dEdPsi/(deltaA2A1.w*lengthPsiCross);
// float4 deltaCross1 = -(deltaCross0+deltaCross2);
f1.xyz += deltaCross0.xyz;
#endif
#ifdef INCLUDE_CHI
// Apply forces based on chi.
float4 ff;
ff.x = (-dEdChi*r)/cross1.w;
ff.y = (deltaA2A1.x*deltaD1A1.x + deltaA2A1.y*deltaD1A1.y + deltaA2A1.z*deltaD1A1.z)/deltaD1A1.w;
ff.z = (deltaD1D2.x*deltaD1A1.x + deltaD1D2.y*deltaD1A1.y + deltaD1D2.z*deltaD1A1.z)/deltaD1A1.w;
ff.w = (dEdChi*r)/cross2.w;
float4 internalF0 = ff.x*cross1;
float4 internalF3 = ff.w*cross2;
float4 s = ff.y*internalF0 - ff.z*internalF3;
f1.xyz -= s.xyz+internalF3.xyz;
f2.xyz += internalF3.xyz;
#endif
#ifdef USE_CUTOFF
}
#endif
}
}
// Write results
int4 bufferIndices = donorBufferIndices[donorIndex];
unsigned int offset1 = atoms.x+bufferIndices.x*PADDED_NUM_ATOMS;
unsigned int offset2 = atoms.y+bufferIndices.y*PADDED_NUM_ATOMS;
float4 force1 = forceBuffers[offset1];
float4 force2 = forceBuffers[offset2];
force1.xyz += f1.xyz;
force2.xyz += f2.xyz;
forceBuffers[offset1] = force1;
forceBuffers[offset2] = force2;
}
energyBuffer[get_global_id(0)] += energy;
}
platforms/opencl/src/kernels/nonbondedExceptions.cl
View file @ b97ff259
......@@ -2,21 +2,16 @@
* Compute nonbonded exceptions.
*/
__kernel void computeNonbondedExceptions(int numAtoms, int numExceptions, float cutoffSquared, float4 periodicBoxSize, __global float4* forceBuffers, __global float* energyBuffer,
__kernel void computeNonbondedExceptions(__global float4* forceBuffers, __global float* energyBuffer,
__global float4* posq, __global float4* params, __global int4* indices) {
int index = get_global_id(0);
float energy = 0.0f;
while (index < numExceptions) {
while (index < NUM_EXCEPTIONS) {
// Look up the data for this bonds.
int4 atoms = indices[index];
float4 exceptionParams = params[index];
float4 delta = posq[atoms.y]-posq[atoms.x];
#ifdef USE_PERIODIC
delta.x -= floor(delta.x/periodicBoxSize.x+0.5f)*periodicBoxSize.x;
delta.y -= floor(delta.y/periodicBoxSize.y+0.5f)*periodicBoxSize.y;
delta.z -= floor(delta.z/periodicBoxSize.z+0.5f)*periodicBoxSize.z;
#endif
// Compute the force.
......@@ -27,27 +22,16 @@ __kernel void computeNonbondedExceptions(int numAtoms, int numExceptions, float
float sig6 = sig2*sig2*sig2;
float dEdR = exceptionParams.z*(12.0f*sig6-6.0f)*sig6;
float tempEnergy = exceptionParams.z*(sig6-1.0f)*sig6;
#ifdef USE_CUTOFF
dEdR += exceptionParams.x*(invR-2.0f*REACTION_FIELD_K*r2);
tempEnergy += exceptionParams.x*(invR+REACTION_FIELD_K*r2-REACTION_FIELD_C);
#else
dEdR += exceptionParams.x*invR;
tempEnergy += exceptionParams.x*invR;
#endif
dEdR *= invR*invR;
#ifdef USE_CUTOFF
if (r2 > cutoffSquared) {
dEdR = 0.0f;
tempEnergy = 0.0f;
}
#endif
tempEnergy += exceptionParams.x*invR;
energy += tempEnergy;
delta.xyz *= dEdR;
// Record the force on each of the two atoms.
unsigned int offsetA = atoms.x+atoms.z*numAtoms;
unsigned int offsetB = atoms.y+atoms.w*numAtoms;
unsigned int offsetA = atoms.x+atoms.z*NUM_ATOMS;
unsigned int offsetB = atoms.y+atoms.w*NUM_ATOMS;
float4 forceA = forceBuffers[offsetA];
float4 forceB = forceBuffers[offsetB];
forceA.xyz -= delta.xyz;
......
platforms/opencl/tests/TestOpenCLEwald.cpp
View file @ b97ff259
......@@ -78,17 +78,20 @@ void testEwaldPME(bool includeExceptions) {
nonbonded->addParticle(1.0, 1.0,0.0);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(-1.0, 1.0,0.0);
if (includeExceptions) {
// Add some exclusions.
for (int i = 0; i < numParticles-1; i++)
nonbonded->addException(i, i+1, 0.0, 1.0, 0.0);
}
system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
vector<Vec3> positions(numParticles);
#include "nacl_amorph.dat"
if (includeExceptions) {
// Add some exclusions.
for (int i = 0; i < numParticles-1; i++) {
Vec3 delta = positions[i]-positions[i+1];
if (sqrt(delta.dot(delta)) < 0.5*cutoff)
nonbonded->addException(i, i+1, i%2 == 0 ? 0.0 : 0.5, 1.0, 0.0);
}
}
// (1) Check whether the Reference and OpenCL platforms agree when using Ewald Method
......
platforms/opencl/tests/TestOpenCLNonbondedForce.cpp
View file @ b97ff259
......@@ -307,10 +307,8 @@ void testCutoff14() {
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces();
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);
force = ONE_4PI_EPS0*q*q*(1.0/(r*r)-2.0*krf*r);
energy = ONE_4PI_EPS0*q*q*(1.0/r+krf*r*r-crf);
force = ONE_4PI_EPS0*q*q/(r*r);
energy = ONE_4PI_EPS0*q*q/r;
if (i == 3) {
force /= 1.2;
energy /= 1.2;
......
platforms/reference/src/ReferenceKernels.cpp
View file @ b97ff259
......@@ -694,7 +694,7 @@ void ReferenceCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
}
if (periodic||ewald||pme)
if (periodic || ewald || pme)
clj.setPeriodic(periodicBoxSize);
if (ewald)
clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
......@@ -703,8 +703,6 @@ void ReferenceCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
ReferenceBondForce refBondForce;
ReferenceLJCoulomb14 nonbonded14;
if (nonbondedMethod == CutoffNonPeriodic || nonbondedMethod == CutoffPeriodic)
nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
}
......@@ -729,8 +727,6 @@ double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
ReferenceBondForce refBondForce;
ReferenceLJCoulomb14 nonbonded14;
if (nonbondedMethod == CutoffNonPeriodic || nonbondedMethod == CutoffPeriodic)
nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
RealOpenMM* energyArray = new RealOpenMM[num14];
for (int i = 0; i < num14; ++i)
energyArray[i] = 0;
......
platforms/reference/src/SimTKReference/ReferenceLJCoulomb14.cpp
View file @ b97ff259
......@@ -37,7 +37,7 @@
--------------------------------------------------------------------------------------- */
ReferenceLJCoulomb14::ReferenceLJCoulomb14( ) : cutoff(false) {
ReferenceLJCoulomb14::ReferenceLJCoulomb14( ) {
// ---------------------------------------------------------------------------------------
......@@ -63,27 +63,6 @@ ReferenceLJCoulomb14::~ReferenceLJCoulomb14( ){
}
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
@param distance the cutoff distance
@param solventDielectric the dielectric constant of the bulk solvent
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceLJCoulomb14::setUseCutoff( RealOpenMM distance, RealOpenMM solventDielectric ) {
cutoff = true;
cutoffDistance = distance;
krf = pow(cutoffDistance, -3.0)*(solventDielectric-1.0)/(2.0*solventDielectric+1.0);
crf = (1.0f/cutoffDistance)*(3.0*solventDielectric)/(2.0*solventDielectric+1.0);
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate LJ 1-4 ixn
......@@ -147,8 +126,6 @@ int ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomC
int atomBIndex = atomIndices[1];
ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0] );
if (cutoff && deltaR[0][ReferenceForce::RIndex] > cutoffDistance)
return ReferenceForce::DefaultReturn;
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM sig2 = inverseR*parameters[0];
......@@ -156,9 +133,6 @@ int ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomC
RealOpenMM sig6 = sig2*sig2*sig2;
RealOpenMM dEdR = parameters[1]*( twelve*sig6 - six )*sig6;
if (cutoff)
dEdR += (RealOpenMM) (ONE_4PI_EPS0*parameters[2]*(inverseR-2.0f*krf*r2));
else
dEdR += (RealOpenMM) (ONE_4PI_EPS0*parameters[2]*inverseR);
dEdR *= inverseR*inverseR;
......@@ -171,9 +145,6 @@ int ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomC
}
RealOpenMM energy = parameters[1]*( sig6 - one )*sig6;
if (cutoff)
energy += (RealOpenMM) (ONE_4PI_EPS0*parameters[2]*(inverseR+krf*r2-crf));
else
energy += (RealOpenMM) (ONE_4PI_EPS0*parameters[2]*inverseR);
// accumulate energies
......
platforms/reference/src/SimTKReference/ReferenceLJCoulomb14.h
View file @ b97ff259
......@@ -31,12 +31,6 @@
class ReferenceLJCoulomb14 : public ReferenceBondIxn {
private:
bool cutoff;
RealOpenMM cutoffDistance;
RealOpenMM krf, crf;
public:
/**---------------------------------------------------------------------------------------
......@@ -55,19 +49,6 @@ class ReferenceLJCoulomb14 : public ReferenceBondIxn {
~ReferenceLJCoulomb14( );
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
@param distance the cutoff distance
@param solventDielectric the dielectric constant of the bulk solvent
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int setUseCutoff( RealOpenMM distance, RealOpenMM solventDielectric );
/**---------------------------------------------------------------------------------------
Calculate Ryckaert-Bellemans bond ixn
......
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
View file @ b97ff259
......@@ -417,10 +417,8 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
int jj = exclusions[i][j];
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
if (r >= cutoffDistance)
continue;
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r;
RealOpenMM dEdR = (RealOpenMM) (ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR);
......
platforms/reference/tests/TestReferenceNonbondedForce.cpp
View file @ b97ff259
......@@ -302,10 +302,8 @@ void testCutoff14() {
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces2 = state.getForces();
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);
force = ONE_4PI_EPS0*q*q*(1.0/(r*r)-2.0*krf*r);
energy = ONE_4PI_EPS0*q*q*(1.0/r+krf*r*r-crf);
force = ONE_4PI_EPS0*q*q/(r*r);
energy = ONE_4PI_EPS0*q*q/r;
if (i == 3) {
force /= 1.2;
energy /= 1.2;
......
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