nonbonded.cu 23.6 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
120
    // used shared memory if the device cannot shuffle
#ifndef ENABLE_SHUFFLE
    __shared__ AtomData localData[THREAD_BLOCK_SIZE];
#endif
    // First loop: process tiles that contain exclusions.
121
122
123
124
125
126
    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;
127
        real3 force = make_real3(0);
128
129
130
131
132
        unsigned int atom1 = x*TILE_SIZE + tgx;
        real4 posq1 = posq[atom1];
        LOAD_ATOM1_PARAMETERS
#ifdef USE_EXCLUSIONS
        tileflags excl = exclusions[pos*TILE_SIZE+tgx];
133
#endif
134
135
136
        const bool hasExclusions = true;
        if (x == y) {
            // This tile is on the diagonal.
137
138
#ifdef ENABLE_SHUFFLE
            real4 shflPosq = posq1;
139
#else
140
141
142
143
144
145
            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
146

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

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

    // Second loop: tiles without exclusions, either from the neighbor list (with cutoff) or just enumerating all
    // of them (no cutoff).
#ifdef USE_CUTOFF
    const unsigned int numTiles = interactionCount[0];
314
315
    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);
316
317
#else
    const unsigned int numTiles = numTileIndices;
318
319
    int pos = (int) (startTileIndex+warp*(long long)numTiles/totalWarps);
    int end = (int) (startTileIndex+(warp+1)*(long long)numTiles/totalWarps);
320
321
322
#endif
    int skipBase = 0;
    int currentSkipIndex = tbx;
323
324
    // atomIndices can probably be shuffled as well
    // but it probably wouldn't make things any faster
325
    __shared__ int atomIndices[THREAD_BLOCK_SIZE];
326
    __shared__ volatile int skipTiles[THREAD_BLOCK_SIZE];
327
328
329
330
331
332
333
334
    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.
335
        int x, y;
336
337
338
        bool singlePeriodicCopy = false;
#ifdef USE_CUTOFF
        if (numTiles <= maxTiles) {
339
            x = tiles[pos];
340
            real4 blockSizeX = blockSize[x];
341
342
343
            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);
344
345
        }
        else
346
#endif
347
        {
348
            y = (int) floor(NUM_BLOCKS+0.5f-SQRT((NUM_BLOCKS+0.5f)*(NUM_BLOCKS+0.5f)-2*pos));
349
350
351
352
353
354
355
356
357
358
359
360
            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;
361
                }
362
363
364
365
                else
                    skipTiles[threadIdx.x] = end;
                skipBase += TILE_SIZE;            
                currentSkipIndex = tbx;
366
            }
367
368
369
370
371
372
            while (skipTiles[currentSkipIndex] < pos)
                currentSkipIndex++;
            includeTile = (skipTiles[currentSkipIndex] != pos);
        }
        if (includeTile) {
            unsigned int atom1 = x*TILE_SIZE + tgx;
373
            // Load atom data for this tile.
374
375
            real4 posq1 = posq[atom1];
            LOAD_ATOM1_PARAMETERS
376
            //const unsigned int localAtomIndex = threadIdx.x;
377
378
379
380
381
382
#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;
383
#ifdef ENABLE_SHUFFLE
384
            DECLARE_LOCAL_PARAMETERS
385
386
387
388
389
390
            real4 shflPosq;
            real3 shflForce;
            shflForce.x = 0.0f;
            shflForce.y = 0.0f;
            shflForce.z = 0.0f;
#endif
391
            if (j < PADDED_NUM_ATOMS) {
392
                // Load position of atom j from from global memory
393
394
#ifdef ENABLE_SHUFFLE
                shflPosq = posq[j];
395
#else
396
397
398
399
400
401
402
403
                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                
404
                LOAD_LOCAL_PARAMETERS_FROM_GLOBAL
405
            }
406
407
408
409
410
411
412
413
414
            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
            }
415
416
417
418
419
#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];
420
                APPLY_PERIODIC_TO_POS_WITH_CENTER(posq1, blockCenterX)
421
#ifdef ENABLE_SHUFFLE
422
                APPLY_PERIODIC_TO_POS_WITH_CENTER(shflPosq, blockCenterX)
423
#else
424
                APPLY_PERIODIC_TO_POS_WITH_CENTER(localData[threadIdx.x], blockCenterX)
425
#endif
426
427
428
                unsigned int tj = tgx;
                for (j = 0; j < TILE_SIZE; j++) {
                    int atom2 = tbx+tj;
429
430
#ifdef ENABLE_SHUFFLE
                    real4 posq2 = shflPosq; 
431
#else
432
433
                    real4 posq2 = make_real4(localData[atom2].x, localData[atom2].y, localData[atom2].z, localData[atom2].q);
#endif
434
                    real3 delta = make_real3(posq2.x-posq1.x, posq2.y-posq1.y, posq2.z-posq1.z);
435
                    real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
436
437
438
439
                    real invR = RSQRT(r2);
                    real r = r2*invR;
                    LOAD_ATOM2_PARAMETERS
                    atom2 = atomIndices[tbx+tj];
440
#ifdef USE_SYMMETRIC
441
                    real dEdR = 0.0f;
442
#else
443
444
                    real3 dEdR1 = make_real3(0);
                    real3 dEdR2 = make_real3(0);
445
446
#endif
#ifdef USE_EXCLUSIONS
447
                    bool isExcluded = (atom1 >= NUM_ATOMS || atom2 >= NUM_ATOMS);
448
#endif
449
450
451
                    real tempEnergy = 0.0f;
                    COMPUTE_INTERACTION
                    energy += tempEnergy;
452
#ifdef INCLUDE_FORCES
453
#ifdef USE_SYMMETRIC
454
455
456
457
                    delta *= dEdR;
                    force.x -= delta.x;
                    force.y -= delta.y;
                    force.z -= delta.z;
458
#ifdef ENABLE_SHUFFLE
459
460
461
                    shflForce.x += delta.x;
                    shflForce.y += delta.y;
                    shflForce.z += delta.z;
462

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

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

            // Write results.
561
#ifdef INCLUDE_FORCES
562
563
564
565
566
567
568
569
570
            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) {
571
572
573
574
#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)));
575
#else
576
577
578
579
                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
580
            }
581
#endif
582
583
        }
        pos++;
584
    }
585
#ifdef INCLUDE_ENERGY
586
    energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy;
587
#endif
588
}