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

Continuing CUDA implementation of CustomNonbondedForce

parent fdc7cc07
...@@ -344,6 +344,7 @@ struct cudaGmxSimulation { ...@@ -344,6 +344,7 @@ struct cudaGmxSimulation {
int4* pCustomExceptionID; // Atom indices for custom nonbonded exceptions int4* pCustomExceptionID; // Atom indices for custom nonbonded exceptions
float4* pCustomExceptionParams; // Parameters for custom nonbonded exceptions float4* pCustomExceptionParams; // Parameters for custom nonbonded exceptions
unsigned int customExceptions; // Number of custom nonbonded exceptions unsigned int customExceptions; // Number of custom nonbonded exceptions
unsigned int customParameters; // Number of parameters for custom nonbonded interactions
float2* pEwaldCosSinSum; // Pointer to the cos/sin sums (ewald) float2* pEwaldCosSinSum; // Pointer to the cos/sin sums (ewald)
unsigned int bonds; // Number of bonds unsigned int bonds; // Number of bonds
int4* pBondID; // Bond atom and output buffer IDs int4* pBondID; // Bond atom and output buffer IDs
......
...@@ -571,6 +571,7 @@ void gpuSetCustomNonbondedParameters(gpuContext gpu, const vector<vector<double> ...@@ -571,6 +571,7 @@ void gpuSetCustomNonbondedParameters(gpuContext gpu, const vector<vector<double>
gpu->sim.nonbondedCutoffSqr = cutoffDistance*cutoffDistance; gpu->sim.nonbondedCutoffSqr = cutoffDistance*cutoffDistance;
gpu->sim.customNonbondedMethod = method; gpu->sim.customNonbondedMethod = method;
gpu->sim.customExceptions = exceptionAtom1.size(); gpu->sim.customExceptions = exceptionAtom1.size();
gpu->sim.customParameters = paramNames.size();
setExclusions(gpu, exclusions); setExclusions(gpu, exclusions);
gpu->psCustomParams = new CUDAStream<float4>(gpu->sim.paddedNumberOfAtoms, 1, "CustomParams"); gpu->psCustomParams = new CUDAStream<float4>(gpu->sim.paddedNumberOfAtoms, 1, "CustomParams");
gpu->sim.pCustomParams = gpu->psCustomParams->_pDevData; gpu->sim.pCustomParams = gpu->psCustomParams->_pDevData;
...@@ -623,6 +624,8 @@ void gpuSetCustomNonbondedParameters(gpuContext gpu, const vector<vector<double> ...@@ -623,6 +624,8 @@ void gpuSetCustomNonbondedParameters(gpuContext gpu, const vector<vector<double>
name << paramNames[i] << j; name << paramNames[i] << j;
combiningRuleParams.push_back(name.str()); combiningRuleParams.push_back(name.str());
} }
for (int i = paramNames.size(); i < 4; i++)
combiningRuleParams.push_back("");
} }
for (int i = 0; i < paramNames.size(); i++) for (int i = 0; i < paramNames.size(); i++)
paramExpressions[i] = createExpression<64>(combiningRules[i], Lepton::Parser::parse(combiningRules[i]).optimize().createProgram(), combiningRuleParams); paramExpressions[i] = createExpression<64>(combiningRules[i], Lepton::Parser::parse(combiningRules[i]).optimize().createProgram(), combiningRuleParams);
......
...@@ -106,27 +106,27 @@ __device__ float kEvaluateExpression_kernel(Expression<SIZE>* expression, float* ...@@ -106,27 +106,27 @@ __device__ float kEvaluateExpression_kernel(Expression<SIZE>* expression, float*
case VARIABLE1: case VARIABLE1:
stack[++stackPointer] = vars1.x; stack[++stackPointer] = vars1.x;
break; break;
// case VARIABLE2: case VARIABLE2:
// stack[++stackPointer] = vars1.y; stack[++stackPointer] = vars1.y;
// break; break;
// case VARIABLE3: case VARIABLE3:
// stack[++stackPointer] = vars1.z; stack[++stackPointer] = vars1.z;
// break; break;
// case VARIABLE4: case VARIABLE4:
// stack[++stackPointer] = vars1.w; stack[++stackPointer] = vars1.w;
// break; break;
case VARIABLE5: case VARIABLE5:
stack[++stackPointer] = vars2.x; stack[++stackPointer] = vars2.x;
break; break;
// case VARIABLE6: case VARIABLE6:
// stack[++stackPointer] = vars2.y; stack[++stackPointer] = vars2.y;
// break; break;
// case VARIABLE7: case VARIABLE7:
// stack[++stackPointer] = vars2.z; stack[++stackPointer] = vars2.z;
// break; break;
// case VARIABLE8: case VARIABLE8:
// stack[++stackPointer] = vars2.w; stack[++stackPointer] = vars2.w;
// break; break;
case ADD: case ADD:
{ {
float temp = stack[stackPointer]; float temp = stack[stackPointer];
......
...@@ -60,7 +60,7 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i ...@@ -60,7 +60,7 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
Atom* psA = &sA[tbx]; Atom* psA = &sA[tbx];
unsigned int i = x + tgx; unsigned int i = x + tgx;
apos = cSim.pPosq[i]; apos = cSim.pPosq[i];
float4 params = make_float4(0,0,0,0); float4 params = cSim.pCustomParams[i];
af.x = 0.0f; af.x = 0.0f;
af.y = 0.0f; af.y = 0.0f;
af.z = 0.0f; af.z = 0.0f;
...@@ -71,44 +71,36 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i ...@@ -71,44 +71,36 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
sA[threadIdx.x].y = apos.y; sA[threadIdx.x].y = apos.y;
sA[threadIdx.x].z = apos.z; sA[threadIdx.x].z = apos.z;
sA[threadIdx.x].params = params; sA[threadIdx.x].params = params;
if (!bExclusionFlag)
{
for (unsigned int j = 0; j < GRID; j++)
{
float4 combinedParams;
float dx = psA[j].x - apos.x;
float dy = psA[j].y - apos.y;
float dz = psA[j].z - apos.z;
#ifdef USE_PERIODIC
dx -= floor(dx/cSim.periodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
dy -= floor(dy/cSim.periodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
dz -= floor(dz/cSim.periodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
#endif
float r = sqrt(dx*dx + dy*dy + dz*dz);
float invR = 1.0f/r;
float dEdR = -kEvaluateExpression_kernel(&forceExp, &stack[MAX_STACK_SIZE*threadIdx.x], r, combinedParams, combinedParams)*invR;
#ifdef USE_CUTOFF
if (r > cSim.nonbondedCutoff)
{
dEdR = 0.0f;
}
#endif
dx *= dEdR;
dy *= dEdR;
dz *= dEdR;
af.x -= dx;
af.y -= dy;
af.z -= dz;
}
}
else // bExclusion
{
unsigned int xi = x>>GRIDBITS; unsigned int xi = x>>GRIDBITS;
unsigned int cell = xi+xi*cSim.paddedNumberOfAtoms/GRID-xi*(xi+1)/2; unsigned int cell = xi+xi*cSim.paddedNumberOfAtoms/GRID-xi*(xi+1)/2;
unsigned int excl = cSim.pExclusion[cSim.pExclusionIndex[cell]+tgx]; unsigned int excl = cSim.pExclusion[cSim.pExclusionIndex[cell]+tgx];
for (unsigned int j = 0; j < GRID; j++) for (unsigned int j = 0; j < GRID; j++)
{ {
float4 combinedParams; // Apply the combining rules to the parameters.
float4 combinedParams = make_float4(0, 0, 0, 0);
for (int k = 0; k < cSim.customParameters; k++)
{
float value = kEvaluateExpression_kernel(&combiningRules[k], &stack[MAX_STACK_SIZE*threadIdx.x], 0.0f, params, psA[j].params);
switch (k)
{
case 0:
combinedParams.x = value;
break;
case 1:
combinedParams.y = value;
break;
case 2:
combinedParams.z = value;
break;
case 3:
combinedParams.w = value;
break;
}
}
// Compute the force.
float dx = psA[j].x - apos.x; float dx = psA[j].x - apos.x;
float dy = psA[j].y - apos.y; float dy = psA[j].y - apos.y;
float dz = psA[j].z - apos.z; float dz = psA[j].z - apos.z;
...@@ -136,7 +128,6 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i ...@@ -136,7 +128,6 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
af.z -= dz; af.z -= dz;
excl >>= 1; excl >>= 1;
} }
}
// Write results // Write results
float4 of; float4 of;
...@@ -166,7 +157,7 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i ...@@ -166,7 +157,7 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
sA[threadIdx.x].x = temp.x; sA[threadIdx.x].x = temp.x;
sA[threadIdx.x].y = temp.y; sA[threadIdx.x].y = temp.y;
sA[threadIdx.x].z = temp.z; sA[threadIdx.x].z = temp.z;
sA[threadIdx.x].params = make_float4(0,0,0,0);; sA[threadIdx.x].params = cSim.pCustomParams[j];
} }
sA[threadIdx.x].fx = 0.0f; sA[threadIdx.x].fx = 0.0f;
sA[threadIdx.x].fy = 0.0f; sA[threadIdx.x].fy = 0.0f;
...@@ -186,7 +177,31 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i ...@@ -186,7 +177,31 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
for (unsigned int j = 0; j < GRID; j++) for (unsigned int j = 0; j < GRID; j++)
{ {
float4 combinedParams; // Apply the combining rules to the parameters.
float4 combinedParams = make_float4(0, 0, 0, 0);
for (int k = 0; k < cSim.customParameters; k++)
{
float value = kEvaluateExpression_kernel(&combiningRules[0], &stack[MAX_STACK_SIZE*threadIdx.x], 0.0f, params, psA[tj].params);
switch (k)
{
case 0:
combinedParams.x = value;
break;
case 1:
combinedParams.y = value;
break;
case 2:
combinedParams.z = value;
break;
case 3:
combinedParams.w = value;
break;
}
}
// Compute the force.
float dx = psA[tj].x - apos.x; float dx = psA[tj].x - apos.x;
float dy = psA[tj].y - apos.y; float dy = psA[tj].y - apos.y;
float dz = psA[tj].z - apos.z; float dz = psA[tj].z - apos.z;
...@@ -225,7 +240,31 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i ...@@ -225,7 +240,31 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
{ {
if ((flags&(1<<j)) != 0) if ((flags&(1<<j)) != 0)
{ {
float4 combinedParams; // Apply the combining rules to the parameters.
float4 combinedParams = make_float4(0, 0, 0, 0);
for (int k = 0; k < cSim.customParameters; k++)
{
float value = kEvaluateExpression_kernel(&combiningRules[0], &stack[MAX_STACK_SIZE*threadIdx.x], 0.0f, params, psA[j].params);
switch (k)
{
case 0:
combinedParams.x = value;
break;
case 1:
combinedParams.y = value;
break;
case 2:
combinedParams.z = value;
break;
case 3:
combinedParams.w = value;
break;
}
}
// Compute the force.
float dx = psA[j].x - apos.x; float dx = psA[j].x - apos.x;
float dy = psA[j].y - apos.y; float dy = psA[j].y - apos.y;
float dz = psA[j].z - apos.z; float dz = psA[j].z - apos.z;
...@@ -300,7 +339,31 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i ...@@ -300,7 +339,31 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
excl = (excl >> tgx) | (excl << (GRID - tgx)); excl = (excl >> tgx) | (excl << (GRID - tgx));
for (unsigned int j = 0; j < GRID; j++) for (unsigned int j = 0; j < GRID; j++)
{ {
float4 combinedParams; // Apply the combining rules to the parameters.
float4 combinedParams = make_float4(0, 0, 0, 0);
for (int k = 0; k < cSim.customParameters; k++)
{
float value = kEvaluateExpression_kernel(&combiningRules[0], &stack[MAX_STACK_SIZE*threadIdx.x], 0.0f, params, psA[tj].params);
switch (k)
{
case 0:
combinedParams.x = value;
break;
case 1:
combinedParams.y = value;
break;
case 2:
combinedParams.z = value;
break;
case 3:
combinedParams.w = value;
break;
}
}
// Compute the force.
float dx = psA[tj].x - apos.x; float dx = psA[tj].x - apos.x;
float dy = psA[tj].y - apos.y; float dy = psA[tj].y - apos.y;
float dz = psA[tj].z - apos.z; float dz = psA[tj].z - apos.z;
......
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