nonbonded.cu 23.5 KB
Newer Older
1
2
#define WARPS_PER_GROUP (THREAD_BLOCK_SIZE/TILE_SIZE)

3
#ifndef ENABLE_SHUFFLE
4
5
6
typedef struct {
    real x, y, z;
    real q;
7
    real fx, fy, fz;
8
9
10
11
12
    ATOM_PARAMETER_DATA
#ifndef PARAMETER_SIZE_IS_EVEN
    real padding;
#endif
} AtomData;
13
14
#endif

15
#ifdef ENABLE_SHUFFLE
16
17
18
19
20
21
//support for 64 bit shuffles
static __inline__ __device__ float real_shfl(float var, int srcLane) {
    return __shfl(var, srcLane);
}

static __inline__ __device__ double real_shfl(double var, int srcLane) {
22
23
24
25
    int hi, lo;
    asm volatile("mov.b64 { %0, %1 }, %2;" : "=r"(lo), "=r"(hi) : "d"(var));
    hi = __shfl(hi, srcLane);
    lo = __shfl(lo, srcLane);
26
27
    return __hiloint2double( hi, lo );
}
28
29
30
31
32
33
34
35
36
37

static __inline__ __device__ long long real_shfl(long long var, int srcLane) {
    int hi, lo;
    asm volatile("mov.b64 { %0, %1 }, %2;" : "=r"(lo), "=r"(hi) : "l"(var));
    hi = __shfl(hi, srcLane);
    lo = __shfl(lo, srcLane);
    // unforunately there isn't an __nv_hiloint2long(hi,lo) intrinsic cast
    int2 fuse; fuse.x = lo; fuse.y = hi;
    return *reinterpret_cast<long long*>(&fuse);
}
38
#endif
39
40

/**
41
42
43
44
45
 * Compute nonbonded interactions. The kernel is separated into two parts,
 * tiles with exclusions and tiles without exclusions. It relies heavily on 
 * implicit warp-level synchronization. A tile is defined by two atom blocks 
 * each of warpsize. Each warp computes a range of tiles.
 * 
46
 * Tiles with exclusions compute the entire set of interactions across
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
 * atom blocks, equal to warpsize*warpsize. In order to avoid access conflicts 
 * the forces are computed and accumulated diagonally in the manner shown below
 * where, suppose
 *
 * [a-h] comprise atom block 1, [i-p] comprise atom block 2
 *
 * 1 denotes the first set of calculations within the warp
 * 2 denotes the second set of calculations within the warp
 * ... etc.
 * 
 *        threads
 *     0 1 2 3 4 5 6 7
 *         atom1 
 * L    a b c d e f g h 
 * o  i 1 2 3 4 5 6 7 8
 * c  j 8 1 2 3 4 5 6 7
 * a  k 7 8 1 2 3 4 5 6
 * l  l 6 7 8 1 2 3 4 5
 * D  m 5 6 7 8 1 2 3 4 
 * a  n 4 5 6 7 8 1 2 3
 * t  o 3 4 5 6 7 8 1 2
68
69
70
 * a  p 2 3 4 5 6 7 8 1
 *
 * Tiles without exclusions read off directly from the neighbourlist interactingAtoms
71
 * and follows the same force accumulation method. If more there are more interactingTiles
72
 * than the size of the neighbourlist initially allocated, the neighbourlist is rebuilt
73
74
75
76
77
78
 * and the full tileset is computed. This should happen on the first step, and very rarely 
 * afterwards.
 *
 * On CUDA devices that support the shuffle intrinsic, on diagonal exclusion tiles use
 * __shfl to broadcast. For all other types of tiles __shfl is used to pass around the 
 * forces, positions, and parameters when computing the forces. 
79
80
81
82
83
84
85
86
 *
 * [out]forceBuffers    - forces on each atom to eventually be accumulated
 * [out]energyBuffer    - energyBuffer to eventually be accumulated
 * [in]posq             - x,y,z,charge 
 * [in]exclusions       - 1024-bit flags denoting atom-atom exclusions for each tile
 * [in]exclusionTiles   - x,y denotes the indices of tiles that have an exclusion
 * [in]startTileIndex   - index into first tile to be processed
 * [in]numTileIndices   - number of tiles this context is responsible for processing
87
 * [in]int tiles        - the atom block for each tile
88
89
90
91
92
93
94
95
96
97
98
99
100
 * [in]interactionCount - total number of tiles that have an interaction
 * [in]maxTiles         - stores the size of the neighbourlist in case it needs 
 *                      - to be expanded
 * [in]periodicBoxSize  - size of the Periodic Box, last dimension (w) not used
 * [in]invPeriodicBox   - inverse of the periodicBoxSize, pre-computed for speed
 * [in]blockCenter      - the center of each block in euclidean coordinates
 * [in]blockSize        - size of the each block, radiating from the center
 *                      - x is half the distance of total length
 *                      - y is half the distance of total width
 *                      - z is half the distance of total height
 *                      - w is not used
 * [in]interactingAtoms - a list of interactions within a given tile     
 *
101
102
 */
extern "C" __global__ void computeNonbonded(
103
        unsigned long long* __restrict__ forceBuffers, mixed* __restrict__ energyBuffer, const real4* __restrict__ posq, const tileflags* __restrict__ exclusions,
104
        const ushort2* __restrict__ exclusionTiles, unsigned int startTileIndex, unsigned int numTileIndices
105
#ifdef USE_CUTOFF
106
107
108
        , const int* __restrict__ tiles, const unsigned int* __restrict__ interactionCount,real4 periodicBoxSize, real4 invPeriodicBoxSize, 
        real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ, unsigned int maxTiles, const real4* __restrict__ blockCenter,
        const real4* __restrict__ blockSize, const unsigned int* __restrict__ interactingAtoms
109
110
#endif
        PARAMETER_ARGUMENTS) {
111
    const unsigned int totalWarps = (blockDim.x*gridDim.x)/TILE_SIZE;
112
113
114
    const unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/TILE_SIZE; // global warpIndex
    const unsigned int tgx = threadIdx.x & (TILE_SIZE-1); // index within the warp
    const unsigned int tbx = threadIdx.x - tgx;           // block warpIndex
115
    mixed energy = 0;
116
117
118
119
    // used shared memory if the device cannot shuffle
#ifndef ENABLE_SHUFFLE
    __shared__ AtomData localData[THREAD_BLOCK_SIZE];
#endif
120

121
    // First loop: process tiles that contain exclusions.
122

123
124
125
126
127
128
    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;
129
        real3 force = make_real3(0);
130
131
132
133
134
        unsigned int atom1 = x*TILE_SIZE + tgx;
        real4 posq1 = posq[atom1];
        LOAD_ATOM1_PARAMETERS
#ifdef USE_EXCLUSIONS
        tileflags excl = exclusions[pos*TILE_SIZE+tgx];
135
#endif
136
137
138
        const bool hasExclusions = true;
        if (x == y) {
            // This tile is on the diagonal.
139
140
#ifdef ENABLE_SHUFFLE
            real4 shflPosq = posq1;
141
#else
142
143
144
145
146
147
            localData[threadIdx.x].x = posq1.x;
            localData[threadIdx.x].y = posq1.y;
            localData[threadIdx.x].z = posq1.z;
            localData[threadIdx.x].q = posq1.w;
            LOAD_LOCAL_PARAMETERS_FROM_1
#endif
148

149
150
            // we do not need to fetch parameters from global since this is a symmetric tile
            // instead we can broadcast the values using shuffle
151
152
            for (unsigned int j = 0; j < TILE_SIZE; j++) {
                int atom2 = tbx+j;
153
                real4 posq2;
154
#ifdef ENABLE_SHUFFLE
155
                BROADCAST_WARP_DATA
156
#else   
157
158
                posq2 = make_real4(localData[atom2].x, localData[atom2].y, localData[atom2].z, localData[atom2].q);
#endif
159
160
                real3 delta = make_real3(posq2.x-posq1.x, posq2.y-posq1.y, posq2.z-posq1.z);
#ifdef USE_PERIODIC
161
                APPLY_PERIODIC_TO_DELTA(delta)
162
163
164
#endif
                real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                real invR = RSQRT(r2);
peastman's avatar
peastman committed
165
                real r = r2*invR;
166
167
168
169
170
171
172
173
                LOAD_ATOM2_PARAMETERS
                atom2 = y*TILE_SIZE+j;
#ifdef USE_SYMMETRIC
                real dEdR = 0.0f;
#else
                real3 dEdR1 = make_real3(0);
                real3 dEdR2 = make_real3(0);
#endif
174
#ifdef USE_EXCLUSIONS
175
176
177
178
179
                bool isExcluded = (atom1 >= NUM_ATOMS || atom2 >= NUM_ATOMS || !(excl & 0x1));
#endif
                real tempEnergy = 0.0f;
                COMPUTE_INTERACTION
                energy += 0.5f*tempEnergy;
180
#ifdef INCLUDE_FORCES
181
182
183
184
#ifdef USE_SYMMETRIC
                force.x -= delta.x*dEdR;
                force.y -= delta.y*dEdR;
                force.z -= delta.z*dEdR;
185
#else
186
187
188
                force.x -= dEdR1.x;
                force.y -= dEdR1.y;
                force.z -= dEdR1.z;
189
#endif
190
#endif
191
#ifdef USE_EXCLUSIONS
192
                excl >>= 1;
193
#endif
194
195
            }
        }
196
197
        else {
            // This is an off-diagonal tile.
198
            unsigned int j = y*TILE_SIZE + tgx;
199
200
201
202
203
204
            real4 shflPosq = posq[j];
#ifdef ENABLE_SHUFFLE
            real3 shflForce;
            shflForce.x = 0.0f;
            shflForce.y = 0.0f;
            shflForce.z = 0.0f;
205
#else
206
207
208
209
210
211
212
213
            localData[threadIdx.x].x = shflPosq.x;
            localData[threadIdx.x].y = shflPosq.y;
            localData[threadIdx.x].z = shflPosq.z;
            localData[threadIdx.x].q = shflPosq.w;
            localData[threadIdx.x].fx = 0.0f;
            localData[threadIdx.x].fy = 0.0f;
            localData[threadIdx.x].fz = 0.0f;
#endif
214
            DECLARE_LOCAL_PARAMETERS
215
            LOAD_LOCAL_PARAMETERS_FROM_GLOBAL
216
#ifdef USE_EXCLUSIONS
217
            excl = (excl >> tgx) | (excl << (TILE_SIZE - tgx));
218
#endif
219
220
221
            unsigned int tj = tgx;
            for (j = 0; j < TILE_SIZE; j++) {
                int atom2 = tbx+tj;
222
223
#ifdef ENABLE_SHUFFLE
                real4 posq2 = shflPosq;
224
#else
225
226
                real4 posq2 = make_real4(localData[atom2].x, localData[atom2].y, localData[atom2].z, localData[atom2].q);
#endif
227
                real3 delta = make_real3(posq2.x-posq1.x, posq2.y-posq1.y, posq2.z-posq1.z);
228
#ifdef USE_PERIODIC
229
                APPLY_PERIODIC_TO_DELTA(delta)
230
231
#endif
                real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
232
233
234
235
                real invR = RSQRT(r2);
                real r = r2*invR;
                LOAD_ATOM2_PARAMETERS
                atom2 = y*TILE_SIZE+tj;
236
#ifdef USE_SYMMETRIC
237
                real dEdR = 0.0f;
238
#else
239
240
                real3 dEdR1 = make_real3(0);
                real3 dEdR2 = make_real3(0);
241
242
#endif
#ifdef USE_EXCLUSIONS
243
                bool isExcluded = (atom1 >= NUM_ATOMS || atom2 >= NUM_ATOMS || !(excl & 0x1));
244
#endif
245
246
247
                real tempEnergy = 0.0f;
                COMPUTE_INTERACTION
                energy += tempEnergy;
248
#ifdef INCLUDE_FORCES
249
#ifdef USE_SYMMETRIC
250
251
252
253
                delta *= dEdR;
                force.x -= delta.x;
                force.y -= delta.y;
                force.z -= delta.z;
254
#ifdef ENABLE_SHUFFLE
255
256
257
                shflForce.x += delta.x;
                shflForce.y += delta.y;
                shflForce.z += delta.z;
258

259
#else
260
261
262
                localData[tbx+tj].fx += delta.x;
                localData[tbx+tj].fy += delta.y;
                localData[tbx+tj].fz += delta.z;
263
264
#endif
#else // !USE_SYMMETRIC
265
266
267
                force.x -= dEdR1.x;
                force.y -= dEdR1.y;
                force.z -= dEdR1.z;
268
#ifdef ENABLE_SHUFFLE
269
270
271
                shflForce.x += dEdR2.x;
                shflForce.y += dEdR2.y;
                shflForce.z += dEdR2.z;
272
#else
273
274
275
                localData[tbx+tj].fx += dEdR2.x;
                localData[tbx+tj].fy += dEdR2.y;
                localData[tbx+tj].fz += dEdR2.z;
276
277
278
279
#endif 
#endif // end USE_SYMMETRIC
#ifdef ENABLE_SHUFFLE
                SHUFFLE_WARP_DATA
280
281
282
283
#endif
#endif
#ifdef USE_EXCLUSIONS
                excl >>= 1;
284
#endif
285
286
                // cycles the indices
                // 0 1 2 3 4 5 6 7 -> 1 2 3 4 5 6 7 0
287
288
                tj = (tj + 1) & (TILE_SIZE - 1);
            }
289
290
            const unsigned int offset = y*TILE_SIZE + tgx;
            // write results for off diagonal tiles
291
#ifdef INCLUDE_FORCES
292
293
294
295
#ifdef ENABLE_SHUFFLE
            atomicAdd(&forceBuffers[offset], static_cast<unsigned long long>((long long) (shflForce.x*0x100000000)));
            atomicAdd(&forceBuffers[offset+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (shflForce.y*0x100000000)));
            atomicAdd(&forceBuffers[offset+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (shflForce.z*0x100000000)));
296
#else
297
298
299
            atomicAdd(&forceBuffers[offset], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fx*0x100000000)));
            atomicAdd(&forceBuffers[offset+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fy*0x100000000)));
            atomicAdd(&forceBuffers[offset+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fz*0x100000000)));
300
#endif
301
#endif
302
        }
303
        // Write results for on and off diagonal tiles
304
#ifdef INCLUDE_FORCES
305
        const unsigned int offset = x*TILE_SIZE + tgx;
306
307
308
        atomicAdd(&forceBuffers[offset], static_cast<unsigned long long>((long long) (force.x*0x100000000)));
        atomicAdd(&forceBuffers[offset+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.y*0x100000000)));
        atomicAdd(&forceBuffers[offset+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.z*0x100000000)));
309
#endif
310
311
312
313
    }

    // Second loop: tiles without exclusions, either from the neighbor list (with cutoff) or just enumerating all
    // of them (no cutoff).
314

315
316
#ifdef USE_CUTOFF
    const unsigned int numTiles = interactionCount[0];
317
318
    if (numTiles > maxTiles)
        return; // There wasn't enough memory for the neighbor list.
319
320
    int pos = (int) (numTiles > maxTiles ? startTileIndex+warp*(long long)numTileIndices/totalWarps : warp*(long long)numTiles/totalWarps);
    int end = (int) (numTiles > maxTiles ? startTileIndex+(warp+1)*(long long)numTileIndices/totalWarps : (warp+1)*(long long)numTiles/totalWarps);
321
322
#else
    const unsigned int numTiles = numTileIndices;
323
324
    int pos = (int) (startTileIndex+warp*(long long)numTiles/totalWarps);
    int end = (int) (startTileIndex+(warp+1)*(long long)numTiles/totalWarps);
325
326
327
#endif
    int skipBase = 0;
    int currentSkipIndex = tbx;
328
329
    // atomIndices can probably be shuffled as well
    // but it probably wouldn't make things any faster
330
    __shared__ int atomIndices[THREAD_BLOCK_SIZE];
331
    __shared__ volatile int skipTiles[THREAD_BLOCK_SIZE];
332
333
334
335
336
337
338
339
    skipTiles[threadIdx.x] = -1;
    
    while (pos < end) {
        const bool hasExclusions = false;
        real3 force = make_real3(0);
        bool includeTile = true;

        // Extract the coordinates of this tile.
340
        int x, y;
341
342
        bool singlePeriodicCopy = false;
#ifdef USE_CUTOFF
343
344
345
346
347
348
349
350
351
352
        x = tiles[pos];
        real4 blockSizeX = blockSize[x];
        singlePeriodicCopy = (0.5f*periodicBoxSize.x-blockSizeX.x >= MAX_CUTOFF &&
                              0.5f*periodicBoxSize.y-blockSizeX.y >= MAX_CUTOFF &&
                              0.5f*periodicBoxSize.z-blockSizeX.z >= MAX_CUTOFF);
#else
        y = (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);
353
            x = (pos-y*NUM_BLOCKS+y*(y+1)/2);
354
        }
355

356
        // Skip over tiles that have exclusions, since they were already processed.
357

358
359
360
361
        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;
362
            }
363
364
365
366
            else
                skipTiles[threadIdx.x] = end;
            skipBase += TILE_SIZE;            
            currentSkipIndex = tbx;
367
        }
368
369
370
371
        while (skipTiles[currentSkipIndex] < pos)
            currentSkipIndex++;
        includeTile = (skipTiles[currentSkipIndex] != pos);
#endif
372
373
        if (includeTile) {
            unsigned int atom1 = x*TILE_SIZE + tgx;
374
            // Load atom data for this tile.
375
376
            real4 posq1 = posq[atom1];
            LOAD_ATOM1_PARAMETERS
377
            //const unsigned int localAtomIndex = threadIdx.x;
378
379
380
381
382
383
#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;
384
#ifdef ENABLE_SHUFFLE
385
            DECLARE_LOCAL_PARAMETERS
386
387
388
389
390
391
            real4 shflPosq;
            real3 shflForce;
            shflForce.x = 0.0f;
            shflForce.y = 0.0f;
            shflForce.z = 0.0f;
#endif
392
            if (j < PADDED_NUM_ATOMS) {
393
                // Load position of atom j from from global memory
394
395
#ifdef ENABLE_SHUFFLE
                shflPosq = posq[j];
396
#else
397
398
399
400
401
402
403
404
                localData[threadIdx.x].x = posq[j].x;
                localData[threadIdx.x].y = posq[j].y;
                localData[threadIdx.x].z = posq[j].z;
                localData[threadIdx.x].q = posq[j].w;
                localData[threadIdx.x].fx = 0.0f;
                localData[threadIdx.x].fy = 0.0f;
                localData[threadIdx.x].fz = 0.0f;
#endif                
405
                LOAD_LOCAL_PARAMETERS_FROM_GLOBAL
406
            }
407
408
409
410
411
412
413
414
415
            else {
#ifdef ENABLE_SHUFFLE
                shflPosq = make_real4(0, 0, 0, 0);
#else
                localData[threadIdx.x].x = 0;
                localData[threadIdx.x].y = 0;
                localData[threadIdx.x].z = 0;
#endif
            }
416
417
418
419
420
#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];
421
                APPLY_PERIODIC_TO_POS_WITH_CENTER(posq1, blockCenterX)
422
#ifdef ENABLE_SHUFFLE
423
                APPLY_PERIODIC_TO_POS_WITH_CENTER(shflPosq, blockCenterX)
424
#else
425
                APPLY_PERIODIC_TO_POS_WITH_CENTER(localData[threadIdx.x], blockCenterX)
426
#endif
427
428
429
                unsigned int tj = tgx;
                for (j = 0; j < TILE_SIZE; j++) {
                    int atom2 = tbx+tj;
430
431
#ifdef ENABLE_SHUFFLE
                    real4 posq2 = shflPosq; 
432
#else
433
434
                    real4 posq2 = make_real4(localData[atom2].x, localData[atom2].y, localData[atom2].z, localData[atom2].q);
#endif
435
                    real3 delta = make_real3(posq2.x-posq1.x, posq2.y-posq1.y, posq2.z-posq1.z);
436
                    real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
437
438
439
440
                    real invR = RSQRT(r2);
                    real r = r2*invR;
                    LOAD_ATOM2_PARAMETERS
                    atom2 = atomIndices[tbx+tj];
441
#ifdef USE_SYMMETRIC
442
                    real dEdR = 0.0f;
443
#else
444
445
                    real3 dEdR1 = make_real3(0);
                    real3 dEdR2 = make_real3(0);
446
447
#endif
#ifdef USE_EXCLUSIONS
448
                    bool isExcluded = (atom1 >= NUM_ATOMS || atom2 >= NUM_ATOMS);
449
#endif
450
451
452
                    real tempEnergy = 0.0f;
                    COMPUTE_INTERACTION
                    energy += tempEnergy;
453
#ifdef INCLUDE_FORCES
454
#ifdef USE_SYMMETRIC
455
456
457
458
                    delta *= dEdR;
                    force.x -= delta.x;
                    force.y -= delta.y;
                    force.z -= delta.z;
459
#ifdef ENABLE_SHUFFLE
460
461
462
                    shflForce.x += delta.x;
                    shflForce.y += delta.y;
                    shflForce.z += delta.z;
463

464
#else
465
466
467
                    localData[tbx+tj].fx += delta.x;
                    localData[tbx+tj].fy += delta.y;
                    localData[tbx+tj].fz += delta.z;
468
469
#endif
#else // !USE_SYMMETRIC
470
471
472
                    force.x -= dEdR1.x;
                    force.y -= dEdR1.y;
                    force.z -= dEdR1.z;
473
#ifdef ENABLE_SHUFFLE
474
475
476
                    shflForce.x += dEdR2.x;
                    shflForce.y += dEdR2.y;
                    shflForce.z += dEdR2.z;
477
#else
478
479
480
                    localData[tbx+tj].fx += dEdR2.x;
                    localData[tbx+tj].fy += dEdR2.y;
                    localData[tbx+tj].fz += dEdR2.z;
481
482
483
#endif 
#endif // end USE_SYMMETRIC
#ifdef ENABLE_SHUFFLE
484
                    SHUFFLE_WARP_DATA
485
#endif
486
#endif
487
                    tj = (tj + 1) & (TILE_SIZE - 1);
488
                }
489
490
            }
            else
491
#endif
492
493
494
495
496
            {
                // We need to apply periodic boundary conditions separately for each interaction.
                unsigned int tj = tgx;
                for (j = 0; j < TILE_SIZE; j++) {
                    int atom2 = tbx+tj;
497
498
#ifdef ENABLE_SHUFFLE
                    real4 posq2 = shflPosq;
499
#else
500
501
                    real4 posq2 = make_real4(localData[atom2].x, localData[atom2].y, localData[atom2].z, localData[atom2].q);
#endif
502
                    real3 delta = make_real3(posq2.x-posq1.x, posq2.y-posq1.y, posq2.z-posq1.z);
503
#ifdef USE_PERIODIC
504
                    APPLY_PERIODIC_TO_DELTA(delta)
505
#endif
506
                    real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
507
508
509
510
                    real invR = RSQRT(r2);
                    real r = r2*invR;
                    LOAD_ATOM2_PARAMETERS
                    atom2 = atomIndices[tbx+tj];
511
#ifdef USE_SYMMETRIC
512
                    real dEdR = 0.0f;
513
#else
514
515
                    real3 dEdR1 = make_real3(0);
                    real3 dEdR2 = make_real3(0);
516
517
#endif
#ifdef USE_EXCLUSIONS
518
                    bool isExcluded = (atom1 >= NUM_ATOMS || atom2 >= NUM_ATOMS);
519
#endif
520
521
522
                    real tempEnergy = 0.0f;
                    COMPUTE_INTERACTION
                    energy += tempEnergy;
523
#ifdef INCLUDE_FORCES
524
#ifdef USE_SYMMETRIC
525
526
527
528
                    delta *= dEdR;
                    force.x -= delta.x;
                    force.y -= delta.y;
                    force.z -= delta.z;
529
#ifdef ENABLE_SHUFFLE
530
531
532
                    shflForce.x += delta.x;
                    shflForce.y += delta.y;
                    shflForce.z += delta.z;
533

534
#else
535
536
537
                    localData[tbx+tj].fx += delta.x;
                    localData[tbx+tj].fy += delta.y;
                    localData[tbx+tj].fz += delta.z;
538
#endif
539
#else // !USE_SYMMETRIC
540
541
542
                    force.x -= dEdR1.x;
                    force.y -= dEdR1.y;
                    force.z -= dEdR1.z;
543
#ifdef ENABLE_SHUFFLE
544
545
546
                    shflForce.x += dEdR2.x;
                    shflForce.y += dEdR2.y;
                    shflForce.z += dEdR2.z;
547
#else
548
549
550
                    localData[tbx+tj].fx += dEdR2.x;
                    localData[tbx+tj].fy += dEdR2.y;
                    localData[tbx+tj].fz += dEdR2.z;
551
552
553
#endif 
#endif // end USE_SYMMETRIC
#ifdef ENABLE_SHUFFLE
554
                    SHUFFLE_WARP_DATA
555
#endif
556
#endif
557
                    tj = (tj + 1) & (TILE_SIZE - 1);
558
559
                }
            }
560
561

            // Write results.
562
#ifdef INCLUDE_FORCES
563
564
565
566
567
568
569
570
571
            atomicAdd(&forceBuffers[atom1], static_cast<unsigned long long>((long long) (force.x*0x100000000)));
            atomicAdd(&forceBuffers[atom1+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.y*0x100000000)));
            atomicAdd(&forceBuffers[atom1+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.z*0x100000000)));
#ifdef USE_CUTOFF
            unsigned int atom2 = atomIndices[threadIdx.x];
#else
            unsigned int atom2 = y*TILE_SIZE + tgx;
#endif
            if (atom2 < PADDED_NUM_ATOMS) {
572
573
574
575
#ifdef ENABLE_SHUFFLE
                atomicAdd(&forceBuffers[atom2], static_cast<unsigned long long>((long long) (shflForce.x*0x100000000)));
                atomicAdd(&forceBuffers[atom2+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (shflForce.y*0x100000000)));
                atomicAdd(&forceBuffers[atom2+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (shflForce.z*0x100000000)));
576
#else
577
578
579
580
                atomicAdd(&forceBuffers[atom2], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fx*0x100000000)));
                atomicAdd(&forceBuffers[atom2+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fy*0x100000000)));
                atomicAdd(&forceBuffers[atom2+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fz*0x100000000)));
#endif
581
            }
582
#endif
583
584
        }
        pos++;
585
    }
586
#ifdef INCLUDE_ENERGY
587
    energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy;
588
#endif
589
}