#define DIELECTRIC_OFFSET 0.009f #define PROBE_RADIUS 0.14f #define SURFACE_AREA_FACTOR -170.351730667551f //-6.0f*3.14159265358979323846f*0.0216f*1000.0f*0.4184f; #define WARPS_PER_GROUP (FORCE_WORK_GROUP_SIZE/TILE_SIZE) /** * Reduce the Born sums to compute the Born radii. */ extern "C" __global__ void reduceBornSum(float alpha, float beta, float gamma, const long long* __restrict__ bornSum, const float2* __restrict__ params, real* __restrict__ bornRadii, real* __restrict__ obcChain) { for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_ATOMS; index += blockDim.x*gridDim.x) { // Get summed Born data real sum = RECIP(0x100000000)*bornSum[index]; // Now calculate Born radius and OBC term. float offsetRadius = params[index].x; sum *= 0.5f*offsetRadius; real sum2 = sum*sum; real sum3 = sum*sum2; real tanhSum = tanh(alpha*sum - beta*sum2 + gamma*sum3); real nonOffsetRadius = offsetRadius + DIELECTRIC_OFFSET; real radius = RECIP(RECIP(offsetRadius) - tanhSum/nonOffsetRadius); real chain = offsetRadius*(alpha - 2.0f*beta*sum + 3.0f*gamma*sum2); chain = (1-tanhSum*tanhSum)*chain / nonOffsetRadius; bornRadii[index] = radius; obcChain[index] = chain; } } /** * Reduce the Born force. */ extern "C" __global__ void reduceBornForce(long long* __restrict__ bornForce, real* __restrict__ energyBuffer, const float2* __restrict__ params, const real* __restrict__ bornRadii, const real* __restrict__ obcChain) { real energy = 0; for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_ATOMS; index += blockDim.x*gridDim.x) { // Get summed Born force real force = RECIP(0x100000000)*bornForce[index]; // Now calculate the actual force float offsetRadius = params[index].x; real bornRadius = bornRadii[index]; real r = offsetRadius+DIELECTRIC_OFFSET+PROBE_RADIUS; real ratio6 = POW((offsetRadius+DIELECTRIC_OFFSET)/bornRadius, 6); real saTerm = SURFACE_AREA_FACTOR*r*r*ratio6; force += saTerm/bornRadius; energy += saTerm; force *= bornRadius*bornRadius*obcChain[index]; bornForce[index] = (long long) (force*0x100000000); } energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy/-6; } typedef struct { real x, y, z; real q; float radius, scaledRadius; real bornSum; } AtomData1; /** * Compute the Born sum. */ extern "C" __global__ void computeBornSum(unsigned long long* __restrict__ global_bornSum, const real4* __restrict__ posq, const float2* __restrict__ global_params, #ifdef USE_CUTOFF const ushort2* __restrict__ tiles, const unsigned int* __restrict__ interactionCount, real4 periodicBoxSize, real4 invPeriodicBoxSize, unsigned int maxTiles, const real4* __restrict__ blockCenter, const unsigned int* __restrict__ interactingAtoms, #else unsigned int numTiles, #endif const ushort2* __restrict__ exclusionTiles) { const unsigned int totalWarps = (blockDim.x*gridDim.x)/TILE_SIZE; const unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/TILE_SIZE; const unsigned int tgx = threadIdx.x & (TILE_SIZE-1); const unsigned int tbx = threadIdx.x - tgx; __shared__ AtomData1 localData[FORCE_WORK_GROUP_SIZE]; // First loop: process tiles that contain exclusions. const unsigned int firstExclusionTile = FIRST_EXCLUSION_TILE+warp*(LAST_EXCLUSION_TILE-FIRST_EXCLUSION_TILE)/totalWarps; const unsigned int lastExclusionTile = FIRST_EXCLUSION_TILE+(warp+1)*(LAST_EXCLUSION_TILE-FIRST_EXCLUSION_TILE)/totalWarps; for (int pos = firstExclusionTile; pos < lastExclusionTile; pos++) { const ushort2 tileIndices = exclusionTiles[pos]; const unsigned int x = tileIndices.x; const unsigned int y = tileIndices.y; real bornSum = 0; unsigned int atom1 = x*TILE_SIZE + tgx; real4 posq1 = posq[atom1]; float2 params1 = global_params[atom1]; if (x == y) { // This tile is on the diagonal. localData[threadIdx.x].x = posq1.x; localData[threadIdx.x].y = posq1.y; localData[threadIdx.x].z = posq1.z; localData[threadIdx.x].q = posq1.w; localData[threadIdx.x].radius = params1.x; localData[threadIdx.x].scaledRadius = params1.y; for (unsigned int j = 0; j < TILE_SIZE; j++) { real3 delta = make_real3(localData[tbx+j].x-posq1.x, localData[tbx+j].y-posq1.y, localData[tbx+j].z-posq1.z); #ifdef USE_PERIODIC delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x; delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y; delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z; #endif real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; #ifdef USE_CUTOFF if (atom1 < NUM_ATOMS && y*TILE_SIZE+j < NUM_ATOMS && r2 < CUTOFF_SQUARED) { #else if (atom1 < NUM_ATOMS && y*TILE_SIZE+j < NUM_ATOMS) { #endif real invR = RSQRT(r2); real r = RECIP(invR); float2 params2 = make_float2(localData[tbx+j].radius, localData[tbx+j].scaledRadius); real rScaledRadiusJ = r+params2.y; if ((j != tgx) && (params1.x < rScaledRadiusJ)) { real l_ij = RECIP(max(params1.x, fabs(r-params2.y))); real u_ij = RECIP(rScaledRadiusJ); real l_ij2 = l_ij*l_ij; real u_ij2 = u_ij*u_ij; real ratio = LOG(u_ij * RECIP(l_ij)); bornSum += l_ij - u_ij + (0.50f*invR*ratio) + 0.25f*(r*(u_ij2-l_ij2) + (params2.y*params2.y*invR)*(l_ij2-u_ij2)); bornSum += (params1.x < params2.y-r ? 2.0f*(RECIP(params1.x)-l_ij) : 0); } } } } else { // This is an off-diagonal tile. unsigned int j = y*TILE_SIZE + tgx; real4 tempPosq = posq[j]; localData[threadIdx.x].x = tempPosq.x; localData[threadIdx.x].y = tempPosq.y; localData[threadIdx.x].z = tempPosq.z; localData[threadIdx.x].q = tempPosq.w; float2 tempParams = global_params[j]; localData[threadIdx.x].radius = tempParams.x; localData[threadIdx.x].scaledRadius = tempParams.y; localData[threadIdx.x].bornSum = 0.0f; // Compute the full set of interactions in this tile. unsigned int tj = tgx; for (j = 0; j < TILE_SIZE; j++) { real3 delta = make_real3(localData[tbx+tj].x-posq1.x, localData[tbx+tj].y-posq1.y, localData[tbx+tj].z-posq1.z); #ifdef USE_PERIODIC delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x; delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y; delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z; #endif real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; #ifdef USE_CUTOFF if (atom1 < NUM_ATOMS && y*TILE_SIZE+tj < NUM_ATOMS && r2 < CUTOFF_SQUARED) { #else if (atom1 < NUM_ATOMS && y*TILE_SIZE+tj < NUM_ATOMS) { #endif real invR = RSQRT(r2); real r = RECIP(invR); float2 params2 = make_float2(localData[tbx+tj].radius, localData[tbx+tj].scaledRadius); real rScaledRadiusJ = r+params2.y; if (params1.x < rScaledRadiusJ) { real l_ij = RECIP(max(params1.x, fabs(r-params2.y))); real u_ij = RECIP(rScaledRadiusJ); real l_ij2 = l_ij*l_ij; real u_ij2 = u_ij*u_ij; real ratio = LOG(u_ij * RECIP(l_ij)); bornSum += l_ij - u_ij + (0.50f*invR*ratio) + 0.25f*(r*(u_ij2-l_ij2) + (params2.y*params2.y*invR)*(l_ij2-u_ij2)); bornSum += (params1.x < params2.y-r ? 2.0f*(RECIP(params1.x)-l_ij) : 0); } real rScaledRadiusI = r+params1.y; if (params2.x < rScaledRadiusI) { real l_ij = RECIP(max(params2.x, fabs(r-params1.y))); real u_ij = RECIP(rScaledRadiusI); real l_ij2 = l_ij*l_ij; real u_ij2 = u_ij*u_ij; real ratio = LOG(u_ij * RECIP(l_ij)); real term = l_ij - u_ij + (0.50f*invR*ratio) + 0.25f*(r*(u_ij2-l_ij2) + (params1.y*params1.y*invR)*(l_ij2-u_ij2)); term += (params2.x < params1.y-r ? 2.0f*(RECIP(params2.x)-l_ij) : 0); localData[tbx+tj].bornSum += term; } } tj = (tj + 1) & (TILE_SIZE - 1); } } // Write results. unsigned int offset = x*TILE_SIZE + tgx; atomicAdd(&global_bornSum[offset], static_cast((long long) (bornSum*0x100000000))); if (x != y) { offset = y*TILE_SIZE + tgx; atomicAdd(&global_bornSum[offset], static_cast((long long) (localData[threadIdx.x].bornSum*0x100000000))); } } // Second loop: tiles without exclusions, either from the neighbor list (with cutoff) or just enumerating all // of them (no cutoff). #ifdef USE_CUTOFF unsigned int numTiles = interactionCount[0]; int pos = warp*(numTiles > maxTiles ? NUM_BLOCKS*(NUM_BLOCKS+1)/2 : numTiles)/totalWarps; int end = (warp+1)*(numTiles > maxTiles ? NUM_BLOCKS*(NUM_BLOCKS+1)/2 : numTiles)/totalWarps; #else int pos = warp*numTiles/totalWarps; int end = (warp+1)*numTiles/totalWarps; #endif int skipBase = 0; int currentSkipIndex = tbx; __shared__ int atomIndices[FORCE_WORK_GROUP_SIZE]; __shared__ volatile int skipTiles[FORCE_WORK_GROUP_SIZE]; skipTiles[threadIdx.x] = -1; while (pos < end) { real bornSum = 0; bool includeTile = true; // Extract the coordinates of this tile. unsigned int x, y; bool singlePeriodicCopy = false; #ifdef USE_CUTOFF if (numTiles <= maxTiles) { ushort2 tileIndices = tiles[pos]; x = tileIndices.x; singlePeriodicCopy = tileIndices.y; } else #endif { y = (unsigned int) floor(NUM_BLOCKS+0.5f-SQRT((NUM_BLOCKS+0.5f)*(NUM_BLOCKS+0.5f)-2*pos)); x = (pos-y*NUM_BLOCKS+y*(y+1)/2); if (x < y || x >= NUM_BLOCKS) { // Occasionally happens due to roundoff error. y += (x < y ? -1 : 1); x = (pos-y*NUM_BLOCKS+y*(y+1)/2); } // Skip over tiles that have exclusions, since they were already processed. while (skipTiles[tbx+TILE_SIZE-1] < pos) { if (skipBase+tgx < NUM_TILES_WITH_EXCLUSIONS) { ushort2 tile = exclusionTiles[skipBase+tgx]; skipTiles[threadIdx.x] = tile.x + tile.y*NUM_BLOCKS - tile.y*(tile.y+1)/2; } else skipTiles[threadIdx.x] = end; skipBase += TILE_SIZE; currentSkipIndex = tbx; } while (skipTiles[currentSkipIndex] < pos) currentSkipIndex++; includeTile = (skipTiles[currentSkipIndex] != pos); } if (includeTile) { unsigned int atom1 = x*TILE_SIZE + tgx; // Load atom data for this tile. real4 posq1 = posq[atom1]; float2 params1 = global_params[atom1]; #ifdef USE_CUTOFF unsigned int j = (numTiles <= maxTiles ? interactingAtoms[pos*TILE_SIZE+tgx] : y*TILE_SIZE + tgx); #else unsigned int j = y*TILE_SIZE + tgx; #endif atomIndices[threadIdx.x] = j; if (j < PADDED_NUM_ATOMS) { real4 tempPosq = posq[j]; localData[threadIdx.x].x = tempPosq.x; localData[threadIdx.x].y = tempPosq.y; localData[threadIdx.x].z = tempPosq.z; localData[threadIdx.x].q = tempPosq.w; float2 tempParams = global_params[j]; localData[threadIdx.x].radius = tempParams.x; localData[threadIdx.x].scaledRadius = tempParams.y; localData[threadIdx.x].bornSum = 0.0f; } #ifdef USE_PERIODIC if (singlePeriodicCopy) { // The box is small enough that we can just translate all the atoms into a single periodic // box, then skip having to apply periodic boundary conditions later. real4 blockCenterX = blockCenter[x]; posq1.x -= floor((posq1.x-blockCenterX.x)*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x; posq1.y -= floor((posq1.y-blockCenterX.y)*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y; posq1.z -= floor((posq1.z-blockCenterX.z)*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z; localData[threadIdx.x].x -= floor((localData[threadIdx.x].x-blockCenterX.x)*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x; localData[threadIdx.x].y -= floor((localData[threadIdx.x].y-blockCenterX.y)*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y; localData[threadIdx.x].z -= floor((localData[threadIdx.x].z-blockCenterX.z)*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z; unsigned int tj = tgx; for (j = 0; j < TILE_SIZE; j++) { real3 delta = make_real3(localData[tbx+tj].x-posq1.x, localData[tbx+tj].y-posq1.y, localData[tbx+tj].z-posq1.z); real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; int atom2 = atomIndices[tbx+tj]; if (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS && r2 < CUTOFF_SQUARED) { real invR = RSQRT(r2); real r = RECIP(invR); float2 params2 = make_float2(localData[tbx+tj].radius, localData[tbx+tj].scaledRadius); real rScaledRadiusJ = r+params2.y; if (params1.x < rScaledRadiusJ) { real l_ij = RECIP(max(params1.x, fabs(r-params2.y))); real u_ij = RECIP(rScaledRadiusJ); real l_ij2 = l_ij*l_ij; real u_ij2 = u_ij*u_ij; real ratio = LOG(u_ij * RECIP(l_ij)); bornSum += l_ij - u_ij + (0.50f*invR*ratio) + 0.25f*(r*(u_ij2-l_ij2) + (params2.y*params2.y*invR)*(l_ij2-u_ij2)); bornSum += (params1.x < params2.y-r ? 2.0f*(RECIP(params1.x)-l_ij) : 0); } real rScaledRadiusI = r+params1.y; if (params2.x < rScaledRadiusI) { real l_ij = RECIP(max(params2.x, fabs(r-params1.y))); real u_ij = RECIP(rScaledRadiusI); real l_ij2 = l_ij*l_ij; real u_ij2 = u_ij*u_ij; real ratio = LOG(u_ij * RECIP(l_ij)); real term = l_ij - u_ij + (0.50f*invR*ratio) + 0.25f*(r*(u_ij2-l_ij2) + (params1.y*params1.y*invR)*(l_ij2-u_ij2)); term += (params2.x < params1.y-r ? 2.0f*(RECIP(params2.x)-l_ij) : 0); localData[tbx+tj].bornSum += term; } } tj = (tj + 1) & (TILE_SIZE - 1); } } else #endif { // We need to apply periodic boundary conditions separately for each interaction. unsigned int tj = tgx; for (j = 0; j < TILE_SIZE; j++) { real3 delta = make_real3(localData[tbx+tj].x-posq1.x, localData[tbx+tj].y-posq1.y, localData[tbx+tj].z-posq1.z); #ifdef USE_PERIODIC delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x; delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y; delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z; #endif real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; int atom2 = atomIndices[tbx+tj]; #ifdef USE_CUTOFF if (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS && r2 < CUTOFF_SQUARED) { #else if (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) { #endif real invR = RSQRT(r2); real r = RECIP(invR); float2 params2 = make_float2(localData[tbx+tj].radius, localData[tbx+tj].scaledRadius); real rScaledRadiusJ = r+params2.y; if (params1.x < rScaledRadiusJ) { real l_ij = RECIP(max(params1.x, fabs(r-params2.y))); real u_ij = RECIP(rScaledRadiusJ); real l_ij2 = l_ij*l_ij; real u_ij2 = u_ij*u_ij; real ratio = LOG(u_ij * RECIP(l_ij)); bornSum += l_ij - u_ij + (0.50f*invR*ratio) + 0.25f*(r*(u_ij2-l_ij2) + (params2.y*params2.y*invR)*(l_ij2-u_ij2)); bornSum += (params1.x < params2.y-r ? 2.0f*(RECIP(params1.x)-l_ij) : 0); } real rScaledRadiusI = r+params1.y; if (params2.x < rScaledRadiusI) { real l_ij = RECIP(max(params2.x, fabs(r-params1.y))); real u_ij = RECIP(rScaledRadiusI); real l_ij2 = l_ij*l_ij; real u_ij2 = u_ij*u_ij; real ratio = LOG(u_ij * RECIP(l_ij)); real term = l_ij - u_ij + (0.50f*invR*ratio) + 0.25f*(r*(u_ij2-l_ij2) + (params1.y*params1.y*invR)*(l_ij2-u_ij2)); term += (params2.x < params1.y-r ? 2.0f*(RECIP(params2.x)-l_ij) : 0); localData[tbx+tj].bornSum += term; } } tj = (tj + 1) & (TILE_SIZE - 1); } } // Write results. atomicAdd(&global_bornSum[atom1], static_cast((long long) (bornSum*0x100000000))); #ifdef USE_CUTOFF unsigned int atom2 = atomIndices[threadIdx.x]; #else unsigned int atom2 = y*TILE_SIZE + tgx; #endif if (atom2 < PADDED_NUM_ATOMS) atomicAdd(&global_bornSum[atom2], static_cast((long long) (localData[threadIdx.x].bornSum*0x100000000))); } pos++; } } typedef struct { real x, y, z; real q; real fx, fy, fz, fw; real bornRadius; } AtomData2; /** * First part of computing the GBSA interaction. */ extern "C" __global__ void computeGBSAForce1(unsigned long long* __restrict__ forceBuffers, unsigned long long* __restrict__ global_bornForce, real* __restrict__ energyBuffer, const real4* __restrict__ posq, const real* __restrict__ global_bornRadii, #ifdef USE_CUTOFF const ushort2* __restrict__ tiles, const unsigned int* __restrict__ interactionCount, real4 periodicBoxSize, real4 invPeriodicBoxSize, unsigned int maxTiles, const real4* __restrict__ blockCenter, const unsigned int* __restrict__ interactingAtoms, #else unsigned int numTiles, #endif const ushort2* __restrict__ exclusionTiles) { const unsigned int totalWarps = (blockDim.x*gridDim.x)/TILE_SIZE; const unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/TILE_SIZE; const unsigned int tgx = threadIdx.x & (TILE_SIZE-1); const unsigned int tbx = threadIdx.x - tgx; real energy = 0; __shared__ AtomData2 localData[FORCE_WORK_GROUP_SIZE]; // First loop: process tiles that contain exclusions. const unsigned int firstExclusionTile = FIRST_EXCLUSION_TILE+warp*(LAST_EXCLUSION_TILE-FIRST_EXCLUSION_TILE)/totalWarps; const unsigned int lastExclusionTile = FIRST_EXCLUSION_TILE+(warp+1)*(LAST_EXCLUSION_TILE-FIRST_EXCLUSION_TILE)/totalWarps; for (int pos = firstExclusionTile; pos < lastExclusionTile; pos++) { const ushort2 tileIndices = exclusionTiles[pos]; const unsigned int x = tileIndices.x; const unsigned int y = tileIndices.y; real4 force = make_real4(0); unsigned int atom1 = x*TILE_SIZE + tgx; real4 posq1 = posq[atom1]; real bornRadius1 = global_bornRadii[atom1]; if (x == y) { // This tile is on the diagonal. localData[threadIdx.x].x = posq1.x; localData[threadIdx.x].y = posq1.y; localData[threadIdx.x].z = posq1.z; localData[threadIdx.x].q = posq1.w; localData[threadIdx.x].bornRadius = bornRadius1; for (unsigned int j = 0; j < TILE_SIZE; j++) { if (atom1 < NUM_ATOMS && y*TILE_SIZE+j < NUM_ATOMS) { real4 posq2 = make_real4(localData[tbx+j].x, localData[tbx+j].y, localData[tbx+j].z, localData[tbx+j].q); real3 delta = make_real3(posq2.x-posq1.x, posq2.y-posq1.y, posq2.z-posq1.z); #ifdef USE_PERIODIC delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x; delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y; delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z; #endif real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; #ifdef USE_CUTOFF if (r2 < CUTOFF_SQUARED) { #endif real invR = RSQRT(r2); real r = RECIP(invR); real bornRadius2 = localData[tbx+j].bornRadius; real alpha2_ij = bornRadius1*bornRadius2; real D_ij = r2*RECIP(4.0f*alpha2_ij); real expTerm = EXP(-D_ij); real denominator2 = r2 + alpha2_ij*expTerm; real denominator = SQRT(denominator2); real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator); real Gpol = tempEnergy*RECIP(denominator2); real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij); real dEdR = Gpol*(1.0f - 0.25f*expTerm); force.w += dGpol_dalpha2_ij*bornRadius2; energy += 0.5f*tempEnergy; delta *= dEdR; force.x -= delta.x; force.y -= delta.y; force.z -= delta.z; #ifdef USE_CUTOFF } #endif } } } else { // This is an off-diagonal tile. unsigned int j = y*TILE_SIZE + tgx; real4 tempPosq = posq[j]; localData[threadIdx.x].x = tempPosq.x; localData[threadIdx.x].y = tempPosq.y; localData[threadIdx.x].z = tempPosq.z; localData[threadIdx.x].q = tempPosq.w; localData[threadIdx.x].bornRadius = global_bornRadii[j]; localData[threadIdx.x].fx = 0.0f; localData[threadIdx.x].fy = 0.0f; localData[threadIdx.x].fz = 0.0f; localData[threadIdx.x].fw = 0.0f; unsigned int tj = tgx; for (j = 0; j < TILE_SIZE; j++) { if (atom1 < NUM_ATOMS && y*TILE_SIZE+tj < NUM_ATOMS) { real4 posq2 = make_real4(localData[tbx+tj].x, localData[tbx+tj].y, localData[tbx+tj].z, localData[tbx+tj].q); real3 delta = make_real3(posq2.x-posq1.x, posq2.y-posq1.y, posq2.z-posq1.z); #ifdef USE_PERIODIC delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x; delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y; delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z; #endif real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; #ifdef USE_CUTOFF if (r2 < CUTOFF_SQUARED) { #endif real invR = RSQRT(r2); real r = RECIP(invR); real bornRadius2 = localData[tbx+tj].bornRadius; real alpha2_ij = bornRadius1*bornRadius2; real D_ij = r2*RECIP(4.0f*alpha2_ij); real expTerm = EXP(-D_ij); real denominator2 = r2 + alpha2_ij*expTerm; real denominator = SQRT(denominator2); real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator); real Gpol = tempEnergy*RECIP(denominator2); real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij); real dEdR = Gpol*(1.0f - 0.25f*expTerm); force.w += dGpol_dalpha2_ij*bornRadius2; energy += tempEnergy; delta *= dEdR; force.x -= delta.x; force.y -= delta.y; force.z -= delta.z; localData[tbx+tj].fx += delta.x; localData[tbx+tj].fy += delta.y; localData[tbx+tj].fz += delta.z; localData[tbx+tj].fw += dGpol_dalpha2_ij*bornRadius1; #ifdef USE_CUTOFF } #endif } tj = (tj + 1) & (TILE_SIZE - 1); } } // Write results. unsigned int offset = x*TILE_SIZE + tgx; atomicAdd(&forceBuffers[offset], static_cast((long long) (force.x*0x100000000))); atomicAdd(&forceBuffers[offset+PADDED_NUM_ATOMS], static_cast((long long) (force.y*0x100000000))); atomicAdd(&forceBuffers[offset+2*PADDED_NUM_ATOMS], static_cast((long long) (force.z*0x100000000))); atomicAdd(&global_bornForce[offset], static_cast((long long) (force.w*0x100000000))); if (x != y) { offset = y*TILE_SIZE + tgx; atomicAdd(&forceBuffers[offset], static_cast((long long) (localData[threadIdx.x].fx*0x100000000))); atomicAdd(&forceBuffers[offset+PADDED_NUM_ATOMS], static_cast((long long) (localData[threadIdx.x].fy*0x100000000))); atomicAdd(&forceBuffers[offset+2*PADDED_NUM_ATOMS], static_cast((long long) (localData[threadIdx.x].fz*0x100000000))); atomicAdd(&global_bornForce[offset], static_cast((long long) (localData[threadIdx.x].fw*0x100000000))); } } // Second loop: tiles without exclusions, either from the neighbor list (with cutoff) or just enumerating all // of them (no cutoff). #ifdef USE_CUTOFF unsigned int numTiles = interactionCount[0]; int pos = warp*(numTiles > maxTiles ? NUM_BLOCKS*(NUM_BLOCKS+1)/2 : numTiles)/totalWarps; int end = (warp+1)*(numTiles > maxTiles ? NUM_BLOCKS*(NUM_BLOCKS+1)/2 : numTiles)/totalWarps; #else int pos = warp*numTiles/totalWarps; int end = (warp+1)*numTiles/totalWarps; #endif int skipBase = 0; int currentSkipIndex = tbx; __shared__ int atomIndices[FORCE_WORK_GROUP_SIZE]; __shared__ volatile int skipTiles[FORCE_WORK_GROUP_SIZE]; skipTiles[threadIdx.x] = -1; while (pos < end) { real4 force = make_real4(0); bool includeTile = true; // Extract the coordinates of this tile. unsigned int x, y; bool singlePeriodicCopy = false; #ifdef USE_CUTOFF if (numTiles <= maxTiles) { ushort2 tileIndices = tiles[pos]; x = tileIndices.x; singlePeriodicCopy = tileIndices.y; } else #endif { y = (unsigned int) floor(NUM_BLOCKS+0.5f-SQRT((NUM_BLOCKS+0.5f)*(NUM_BLOCKS+0.5f)-2*pos)); x = (pos-y*NUM_BLOCKS+y*(y+1)/2); if (x < y || x >= NUM_BLOCKS) { // Occasionally happens due to roundoff error. y += (x < y ? -1 : 1); x = (pos-y*NUM_BLOCKS+y*(y+1)/2); } // Skip over tiles that have exclusions, since they were already processed. while (skipTiles[tbx+TILE_SIZE-1] < pos) { if (skipBase+tgx < NUM_TILES_WITH_EXCLUSIONS) { ushort2 tile = exclusionTiles[skipBase+tgx]; skipTiles[threadIdx.x] = tile.x + tile.y*NUM_BLOCKS - tile.y*(tile.y+1)/2; } else skipTiles[threadIdx.x] = end; skipBase += TILE_SIZE; currentSkipIndex = tbx; } while (skipTiles[currentSkipIndex] < pos) currentSkipIndex++; includeTile = (skipTiles[currentSkipIndex] != pos); } if (includeTile) { unsigned int atom1 = x*TILE_SIZE + tgx; // Load atom data for this tile. real4 posq1 = posq[atom1]; real bornRadius1 = global_bornRadii[atom1]; #ifdef USE_CUTOFF unsigned int j = (numTiles <= maxTiles ? interactingAtoms[pos*TILE_SIZE+tgx] : y*TILE_SIZE + tgx); #else unsigned int j = y*TILE_SIZE + tgx; #endif atomIndices[threadIdx.x] = j; if (j < PADDED_NUM_ATOMS) { real4 tempPosq = posq[j]; localData[threadIdx.x].x = tempPosq.x; localData[threadIdx.x].y = tempPosq.y; localData[threadIdx.x].z = tempPosq.z; localData[threadIdx.x].q = tempPosq.w; localData[threadIdx.x].bornRadius = global_bornRadii[j]; localData[threadIdx.x].fx = 0.0f; localData[threadIdx.x].fy = 0.0f; localData[threadIdx.x].fz = 0.0f; localData[threadIdx.x].fw = 0.0f; } #ifdef USE_PERIODIC if (singlePeriodicCopy) { // The box is small enough that we can just translate all the atoms into a single periodic // box, then skip having to apply periodic boundary conditions later. real4 blockCenterX = blockCenter[x]; posq1.x -= floor((posq1.x-blockCenterX.x)*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x; posq1.y -= floor((posq1.y-blockCenterX.y)*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y; posq1.z -= floor((posq1.z-blockCenterX.z)*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z; localData[threadIdx.x].x -= floor((localData[threadIdx.x].x-blockCenterX.x)*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x; localData[threadIdx.x].y -= floor((localData[threadIdx.x].y-blockCenterX.y)*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y; localData[threadIdx.x].z -= floor((localData[threadIdx.x].z-blockCenterX.z)*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z; unsigned int tj = tgx; for (j = 0; j < TILE_SIZE; j++) { int atom2 = atomIndices[tbx+tj]; if (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) { real4 posq2 = make_real4(localData[tbx+tj].x, localData[tbx+tj].y, localData[tbx+tj].z, localData[tbx+tj].q); real3 delta = make_real3(posq2.x-posq1.x, posq2.y-posq1.y, posq2.z-posq1.z); real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; if (r2 < CUTOFF_SQUARED) { real invR = RSQRT(r2); real r = RECIP(invR); real bornRadius2 = localData[tbx+tj].bornRadius; real alpha2_ij = bornRadius1*bornRadius2; real D_ij = r2*RECIP(4.0f*alpha2_ij); real expTerm = EXP(-D_ij); real denominator2 = r2 + alpha2_ij*expTerm; real denominator = SQRT(denominator2); real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator); real Gpol = tempEnergy*RECIP(denominator2); real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij); real dEdR = Gpol*(1.0f - 0.25f*expTerm); force.w += dGpol_dalpha2_ij*bornRadius2; energy += tempEnergy; delta *= dEdR; force.x -= delta.x; force.y -= delta.y; force.z -= delta.z; localData[tbx+tj].fx += delta.x; localData[tbx+tj].fy += delta.y; localData[tbx+tj].fz += delta.z; localData[tbx+tj].fw += dGpol_dalpha2_ij*bornRadius1; } } tj = (tj + 1) & (TILE_SIZE - 1); } } else #endif { // We need to apply periodic boundary conditions separately for each interaction. unsigned int tj = tgx; for (j = 0; j < TILE_SIZE; j++) { int atom2 = atomIndices[tbx+tj]; if (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) { real4 posq2 = make_real4(localData[tbx+tj].x, localData[tbx+tj].y, localData[tbx+tj].z, localData[tbx+tj].q); real3 delta = make_real3(posq2.x-posq1.x, posq2.y-posq1.y, posq2.z-posq1.z); #ifdef USE_PERIODIC delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x; delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y; delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z; #endif real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; #ifdef USE_CUTOFF if (r2 < CUTOFF_SQUARED) { #endif real invR = RSQRT(r2); real r = RECIP(invR); real bornRadius2 = localData[tbx+tj].bornRadius; real alpha2_ij = bornRadius1*bornRadius2; real D_ij = r2*RECIP(4.0f*alpha2_ij); real expTerm = EXP(-D_ij); real denominator2 = r2 + alpha2_ij*expTerm; real denominator = SQRT(denominator2); real tempEnergy = (PREFACTOR*posq1.w*posq2.w)*RECIP(denominator); real Gpol = tempEnergy*RECIP(denominator2); real dGpol_dalpha2_ij = -0.5f*Gpol*expTerm*(1.0f+D_ij); real dEdR = Gpol*(1.0f - 0.25f*expTerm); force.w += dGpol_dalpha2_ij*bornRadius2; energy += tempEnergy; delta *= dEdR; force.x -= delta.x; force.y -= delta.y; force.z -= delta.z; localData[tbx+tj].fx += delta.x; localData[tbx+tj].fy += delta.y; localData[tbx+tj].fz += delta.z; localData[tbx+tj].fw += dGpol_dalpha2_ij*bornRadius1; #ifdef USE_CUTOFF } #endif } tj = (tj + 1) & (TILE_SIZE - 1); } } // Write results. atomicAdd(&forceBuffers[atom1], static_cast((long long) (force.x*0x100000000))); atomicAdd(&forceBuffers[atom1+PADDED_NUM_ATOMS], static_cast((long long) (force.y*0x100000000))); atomicAdd(&forceBuffers[atom1+2*PADDED_NUM_ATOMS], static_cast((long long) (force.z*0x100000000))); atomicAdd(&global_bornForce[atom1], static_cast((long long) (force.w*0x100000000))); #ifdef USE_CUTOFF unsigned int atom2 = atomIndices[threadIdx.x]; #else unsigned int atom2 = y*TILE_SIZE + tgx; #endif if (atom2 < PADDED_NUM_ATOMS) { atomicAdd(&forceBuffers[atom2], static_cast((long long) (localData[threadIdx.x].fx*0x100000000))); atomicAdd(&forceBuffers[atom2+PADDED_NUM_ATOMS], static_cast((long long) (localData[threadIdx.x].fy*0x100000000))); atomicAdd(&forceBuffers[atom2+2*PADDED_NUM_ATOMS], static_cast((long long) (localData[threadIdx.x].fz*0x100000000))); atomicAdd(&global_bornForce[atom2], static_cast((long long) (localData[threadIdx.x].fw*0x100000000))); } } pos++; } energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy; }