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

OpenCL version of CustomHbondForce now supports exclusions

parent 6ba9a156
......@@ -2530,6 +2530,10 @@ OpenCLCalcCustomHbondForceKernel::~OpenCLCalcCustomHbondForceKernel() {
delete acceptorBufferIndices;
if (globals != NULL)
delete globals;
if (donorExclusions != NULL)
delete donorExclusions;
if (acceptorExclusions != NULL)
delete acceptorExclusions;
if (tabulatedFunctionParams != NULL)
delete tabulatedFunctionParams;
for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
......@@ -2607,12 +2611,36 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
// Record exclusions.
vector<vector<int> > exclusionList(numDonors);
vector<mm_int4> donorExclusionVector(numDonors, mm_int4(-1, -1, -1, -1));
vector<mm_int4> acceptorExclusionVector(numAcceptors, mm_int4(-1, -1, -1, -1));
for (int i = 0; i < force.getNumExclusions(); i++) {
int donor, acceptor;
force.getExclusionParticles(i, donor, acceptor);
exclusionList[donor].push_back(acceptor);
if (donorExclusionVector[donor].x == -1)
donorExclusionVector[donor].x = acceptor;
else if (donorExclusionVector[donor].y == -1)
donorExclusionVector[donor].y = acceptor;
else if (donorExclusionVector[donor].z == -1)
donorExclusionVector[donor].z = acceptor;
else if (donorExclusionVector[donor].w == -1)
donorExclusionVector[donor].w = acceptor;
else
throw OpenMMException("CustomHbondForce: OpenCLPlatform does not support more than four exclusions per donor");
if (acceptorExclusionVector[acceptor].x == -1)
acceptorExclusionVector[acceptor].x = donor;
else if (acceptorExclusionVector[acceptor].y == -1)
acceptorExclusionVector[acceptor].y = donor;
else if (acceptorExclusionVector[acceptor].z == -1)
acceptorExclusionVector[acceptor].z = donor;
else if (acceptorExclusionVector[acceptor].w == -1)
acceptorExclusionVector[acceptor].w = donor;
else
throw OpenMMException("CustomHbondForce: OpenCLPlatform does not support more than four exclusions per acceptor");
}
donorExclusions = new OpenCLArray<mm_int4>(cl, numDonors, "customHbondDonorExclusions");
acceptorExclusions = new OpenCLArray<mm_int4>(cl, numDonors, "customHbondAcceptorExclusions");
donorExclusions->upload(donorExclusionVector);
acceptorExclusions->upload(acceptorExclusionVector);
// Record the tabulated functions.
......@@ -2827,7 +2855,8 @@ void OpenCLCalcCustomHbondForceKernel::initialize(const System& system, const Cu
}
if (force.getNonbondedMethod() != CustomHbondForce::NoCutoff && force.getNonbondedMethod() != CustomHbondForce::CutoffNonPeriodic)
defines["USE_PERIODIC"] = "1";
if (force.getNumExclusions() > 0)
defines["USE_EXCLUSIONS"] = "1";
cl::Program program = cl.createProgram(cl.replaceStrings(OpenCLKernelSources::customHbondForce, replacements), defines);
donorKernel = cl::Kernel(program, "computeDonorForces");
acceptorKernel = cl::Kernel(program, "computeAcceptorForces");
......@@ -2851,6 +2880,7 @@ void OpenCLCalcCustomHbondForceKernel::executeForces(ContextImpl& context) {
donorKernel.setArg<cl::Buffer>(index++, cl.getForceBuffers().getDeviceBuffer());
donorKernel.setArg<cl::Buffer>(index++, cl.getEnergyBuffer().getDeviceBuffer());
donorKernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
donorKernel.setArg<cl::Buffer>(index++, donorExclusions->getDeviceBuffer());
donorKernel.setArg<cl::Buffer>(index++, donors->getDeviceBuffer());
donorKernel.setArg<cl::Buffer>(index++, acceptors->getDeviceBuffer());
donorKernel.setArg<cl::Buffer>(index++, donorBufferIndices->getDeviceBuffer());
......@@ -2874,6 +2904,7 @@ void OpenCLCalcCustomHbondForceKernel::executeForces(ContextImpl& context) {
acceptorKernel.setArg<cl::Buffer>(index++, cl.getForceBuffers().getDeviceBuffer());
acceptorKernel.setArg<cl::Buffer>(index++, cl.getEnergyBuffer().getDeviceBuffer());
acceptorKernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
acceptorKernel.setArg<cl::Buffer>(index++, acceptorExclusions->getDeviceBuffer());
acceptorKernel.setArg<cl::Buffer>(index++, donors->getDeviceBuffer());
acceptorKernel.setArg<cl::Buffer>(index++, acceptors->getDeviceBuffer());
acceptorKernel.setArg<cl::Buffer>(index++, acceptorBufferIndices->getDeviceBuffer());
......
......@@ -671,7 +671,8 @@ 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) {
donorBufferIndices(NULL), acceptorBufferIndices(NULL), globals(NULL), donorExclusions(NULL), acceptorExclusions(NULL),
tabulatedFunctionParams(NULL), system(system) {
}
~OpenCLCalcCustomHbondForceKernel();
/**
......@@ -705,6 +706,8 @@ private:
OpenCLArray<mm_int4>* acceptors;
OpenCLArray<mm_int4>* donorBufferIndices;
OpenCLArray<mm_int4>* acceptorBufferIndices;
OpenCLArray<mm_int4>* donorExclusions;
OpenCLArray<mm_int4>* acceptorExclusions;
OpenCLArray<mm_float4>* tabulatedFunctionParams;
std::vector<std::string> globalParamNames;
std::vector<cl_float> globalParamValues;
......
......@@ -55,8 +55,8 @@ float4 computeCross(float4 vec1, float4 vec2) {
/**
* Compute forces on donors.
*/
__kernel void computeDonorForces(__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, __local float4* posBuffer
__kernel void computeDonorForces(__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global int4* exclusions,
__global int4* donorAtoms, __global int4* acceptorAtoms, __global int4* donorBufferIndices, __local float4* posBuffer
PARAMETER_ARGUMENTS) {
float energy = 0.0f;
float4 f1 = 0;
......@@ -66,13 +66,16 @@ __kernel void computeDonorForces(__global float4* forceBuffers, __global float*
// Load information about the donor this thread will compute forces on.
int donorIndex = donorStart+get_global_id(0);
int4 atoms;
int4 atoms, exclusionIndices;
float4 d1, d2, d3;
if (donorIndex < NUM_DONORS) {
atoms = donorAtoms[donorIndex];
d1 = posq[atoms.x];
d2 = posq[atoms.y];
d3 = posq[atoms.z];
#ifdef USE_EXCLUSIONS
exclusionIndices = exclusions[donorIndex];
#endif
}
else
atoms = (int4) (-1, -1, -1, -1);
......@@ -89,6 +92,11 @@ __kernel void computeDonorForces(__global float4* forceBuffers, __global float*
barrier(CLK_LOCAL_MEM_FENCE);
if (donorIndex < NUM_DONORS) {
for (int index = 0; index < blockSize; index++) {
#ifdef USE_EXCLUSIONS
int acceptorIndex = acceptorStart+index;
if (acceptorIndex == exclusionIndices.x || acceptorIndex == exclusionIndices.y || acceptorIndex == exclusionIndices.z || acceptorIndex == exclusionIndices.w)
continue;
#endif
// Compute the interaction between a donor and an acceptor.
float4 a1 = posBuffer[3*index];
......@@ -133,8 +141,8 @@ __kernel void computeDonorForces(__global float4* forceBuffers, __global float*
/**
* Compute forces on acceptors.
*/
__kernel void computeAcceptorForces(__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* acceptorBufferIndices, __local float4* posBuffer
__kernel void computeAcceptorForces(__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global int4* exclusions,
__global int4* donorAtoms, __global int4* acceptorAtoms, __global int4* acceptorBufferIndices, __local float4* posBuffer
PARAMETER_ARGUMENTS) {
float4 f1 = 0;
float4 f2 = 0;
......@@ -143,13 +151,16 @@ __kernel void computeAcceptorForces(__global float4* forceBuffers, __global floa
// Load information about the acceptor this thread will compute forces on.
int acceptorIndex = acceptorStart+get_global_id(0);
int4 atoms;
int4 atoms, exclusionIndices;
float4 a1, a2, a3;
if (acceptorIndex < NUM_ACCEPTORS) {
atoms = acceptorAtoms[acceptorIndex];
a1 = posq[atoms.x];
a2 = posq[atoms.y];
a3 = posq[atoms.z];
#ifdef USE_EXCLUSIONS
exclusionIndices = exclusions[acceptorIndex];
#endif
}
else
atoms = (int4) (-1, -1, -1, -1);
......@@ -166,6 +177,11 @@ __kernel void computeAcceptorForces(__global float4* forceBuffers, __global floa
barrier(CLK_LOCAL_MEM_FENCE);
if (acceptorIndex < NUM_ACCEPTORS) {
for (int index = 0; index < blockSize; index++) {
#ifdef USE_EXCLUSIONS
int donorIndex = donorStart+index;
if (donorIndex == exclusionIndices.x || donorIndex == exclusionIndices.y || donorIndex == exclusionIndices.z || donorIndex == exclusionIndices.w)
continue;
#endif
// Compute the interaction between a donor and an acceptor.
float4 d1 = posBuffer[3*index];
......
......@@ -212,7 +212,7 @@ void testCustomFunctions() {
int main() {
try {
testHbond();
// testExclusions();
testExclusions();
testCutoff();
testCustomFunctions();
}
......
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