kCalculateCDLJForces.h 25 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2009 Stanford University and the Authors.           *
 * Authors: Scott Le Grand, Peter Eastman                                     *
 * Contributors:                                                              *
 *                                                                            *
13
14
15
16
 * This program is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU Lesser General Public License as published   *
 * by the Free Software Foundation, either version 3 of the License, or       *
 * (at your option) any later version.                                        *
17
 *                                                                            *
18
19
20
21
 * This program is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
 * GNU Lesser General Public License for more details.                        *
22
 *                                                                            *
23
24
 * You should have received a copy of the GNU Lesser General Public License   *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
25
26
27
28
29
30
31
32
 * -------------------------------------------------------------------------- */

/**
 * This file contains the kernels for evalauating nonbonded forces.  It is included
 * several times in kCalculateCDLJForces.cu with different #defines to generate
 * different versions of the kernels.
 */

33
__global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUnit, float * gpu_energies)
34
35
36
37
{
    extern __shared__ Atom sA[];
    unsigned int totalWarps = cSim.nonbond_blocks*cSim.nonbond_threads_per_block/GRID;
    unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
38
    unsigned int numWorkUnits = cSim.pInteractionCount[0];
39
40
    unsigned int pos = warp*numWorkUnits/totalWarps;
    unsigned int end = (warp+1)*numWorkUnits/totalWarps;
41
42
43
    unsigned int energyIndex = blockIdx.x * blockDim.x + threadIdx.x;
    float CDLJ_energy;
    float energy = 0.0f;
44
45
46
#ifdef USE_CUTOFF
    float3* tempBuffer = (float3*) &sA[cSim.nonbond_threads_per_block];
#endif
47

48
#ifdef USE_EWALD
49
50
    float PI           = 3.14159265358979323846f;
    float SQRT_PI      = sqrt(PI);
51
52
#endif

53
    unsigned int lasty = 0xFFFFFFFF;
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
    while (pos < end)
    {

        // Extract cell coordinates from appropriate work unit
        unsigned int x = workUnit[pos];
        unsigned int y = ((x >> 2) & 0x7fff) << GRIDBITS;
        bool bExclusionFlag = (x & 0x1);
        x = (x >> 17) << GRIDBITS;
        float4      apos;   // Local atom x, y, z, q
        float3      af;     // Local atom fx, fy, fz
        float dx;
        float dy;
        float dz;
        float r2;
        float invR;
        float sig;
        float sig2;
        float sig6;
        float eps;
        float dEdR;
        unsigned int tgx = threadIdx.x & (GRID - 1);
        unsigned int tbx = threadIdx.x - tgx;
76
        unsigned int tj = tgx;
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
        Atom* psA = &sA[tbx];
        unsigned int i      = x + tgx;
        apos                = cSim.pPosq[i];
        float2 a            = cSim.pAttr[i];
        af.x                = 0.0f;
        af.y                = 0.0f;
        af.z                = 0.0f;
        if (x == y) // Handle diagonals uniquely at 50% efficiency
        {
            // Read fixed atom data into registers and GRF
            sA[threadIdx.x].x   = apos.x;
            sA[threadIdx.x].y   = apos.y;
            sA[threadIdx.x].z   = apos.z;
            sA[threadIdx.x].q   = apos.w;
            sA[threadIdx.x].sig = a.x;
            sA[threadIdx.x].eps = a.y;
            apos.w             *= cSim.epsfac;
            if (!bExclusionFlag)
            {
                for (unsigned int j = 0; j < GRID; j++)
                {
                    dx              = psA[j].x - apos.x;
                    dy              = psA[j].y - apos.y;
                    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
                    r2              = dx * dx + dy * dy + dz * dz;
                    invR            = 1.0f / sqrt(r2);
                    sig             = a.x + psA[j].sig;
                    sig2            = invR * sig;
                    sig2           *= sig2;
                    sig6            = sig2 * sig2 * sig2;
                    eps             = a.y * psA[j].eps;
                    dEdR            = eps * (12.0f * sig6 - 6.0f) * sig6;
114
115
		    /* E */ 
		    CDLJ_energy     = eps * (sig6 - 1.0f) * sig6;
116
#ifdef USE_CUTOFF
117
    #ifdef USE_EWALD
118
                    float r         = sqrt(r2);
119
                    float alphaR    = cSim.alphaEwald * r;
120
121
122
                    dEdR           += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
		    /* E */
		    CDLJ_energy    += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
123
    #else
124
                    dEdR           += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2);
125
126
		    /* E */
		    CDLJ_energy    += apos.w * psA[j].q * (invR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
127
    #endif
128
129
#else
                    dEdR           += apos.w * psA[j].q * invR;
130
131
		    /* E */
		    CDLJ_energy    += apos.w * psA[j].q * invR;
132
133
134
135
136
137
#endif
                    dEdR           *= invR * invR;
#ifdef USE_CUTOFF
                    if (r2 > cSim.nonbondedCutoffSqr)
                    {
                        dEdR = 0.0f;
138
139
			/* E */
			CDLJ_energy = 0.0f;
140
141
                    }
#endif
142
143
		    /* E */
		    energy         += CDLJ_energy / 2.0;
144
145
146
147
148
149
150
151
152
153
                    dx             *= dEdR;
                    dy             *= dEdR;
                    dz             *= dEdR;
                    af.x           -= dx;
                    af.y           -= dy;
                    af.z           -= dz;
                }
            }
            else  // bExclusion
            {
154
                unsigned int xi   = x>>GRIDBITS;
155
                unsigned int cell          = xi+xi*cSim.paddedNumberOfAtoms/GRID-xi*(xi+1)/2;
156
                unsigned int excl = cSim.pExclusion[cSim.pExclusionIndex[cell]+tgx];
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
                for (unsigned int j = 0; j < GRID; j++)
                {
                    dx              = psA[j].x - apos.x;
                    dy              = psA[j].y - apos.y;
                    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
                    r2              = dx * dx + dy * dy + dz * dz;
                    invR            = 1.0f / sqrt(r2);
                    sig             = a.x + psA[j].sig;
                    sig2            = invR * sig;
                    sig2           *= sig2;
                    sig6            = sig2 * sig2 * sig2;
                    eps             = a.y * psA[j].eps;
                    dEdR            = eps * (12.0f * sig6 - 6.0f) * sig6;
175
176
		    /* E */
		    CDLJ_energy     = eps * (sig6 - 1.0f) * sig6;
177
#ifdef USE_CUTOFF
178
    #ifdef USE_EWALD
179
                    float r         = sqrt(r2);
180
                    float alphaR    = cSim.alphaEwald * r;
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
                    dEdR           += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
		    /* E */
		    //CDLJ_energy    += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
                    // direct space sum
                    // CDLJ_energy      += 0.5f * apos.w * psA[j].q * erfc(alphaR) * invR;
                    // reciprocal space sum
                    CDLJ_energy     = 0.5f * apos.w * psA[j].q / PI;
                    CDLJ_energy    *= exp ( - alphaR * alphaR) / SQRT_PI;
                    // self Ewald Energy
                    // factor 2 necessary here, because energy is divided by 2 below
//                    CDLJ_energy     += -cSim.epsfac * (psA[i].q * psA[i].q) * cSim.alphaEwald / SQRT_PI;

                    //    if (r != 0) energy = cSim.alphaEwald;
                   // if (r != 0) energy = apos.w;
   #else
196
                    dEdR           += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2);
197
198
                    /* E */
		    CDLJ_energy    += apos.w * psA[j].q * (invR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
199
    #endif
200
201
#else
                    dEdR           += apos.w * psA[j].q * invR;
202
203
                    /* E */
		    CDLJ_energy    += apos.w * psA[j].q * invR;
204
205
206
207
208
209
210
211
212
#endif
                    dEdR           *= invR * invR;
#ifdef USE_CUTOFF
                    if (!(excl & 0x1) || r2 > cSim.nonbondedCutoffSqr)
#else
                    if (!(excl & 0x1))
#endif
                    {
                        dEdR = 0.0f;
213
214
			/* E */
		        CDLJ_energy  = 0.0f;
215
                    }
216
217
		    /* E */
                    energy         += CDLJ_energy / 2.0;
218
219
220
221
222
223
224
225
226
227
228
229
230
                    dx             *= dEdR;
                    dy             *= dEdR;
                    dz             *= dEdR;
                    af.x           -= dx;
                    af.y           -= dy;
                    af.z           -= dz;
                    excl          >>= 1;
                }
            }

            // Write results
            float4 of;
#ifdef USE_OUTPUT_BUFFER_PER_WARP
231
            unsigned int offset                          = x + tgx + warp*cSim.stride;
232
233
234
235
236
237
238
239
240
241
            of                                  = cSim.pForce4a[offset];
            of.x                               += af.x;
            of.y                               += af.y;
            of.z                               += af.z;
            cSim.pForce4a[offset]               = of;
#else
            of.x                                = af.x;
            of.y                                = af.y;
            of.z                                = af.z;
            of.w                                = 0.0f;
242
            unsigned int offset                          = x + tgx + (x >> GRIDBITS) * cSim.stride;
243
244
245
246
247
248
249
250
            cSim.pForce4a[offset]               = of;
#endif
        }
        else        // 100% utilization
        {
            // Read fixed atom data into registers and GRF
            if (lasty != y)
            {
251
                unsigned int j                   = y + tgx;
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
                float4 temp             = cSim.pPosq[j];
                float2 temp1            = cSim.pAttr[j];
                sA[threadIdx.x].x       = temp.x;
                sA[threadIdx.x].y       = temp.y;
                sA[threadIdx.x].z       = temp.z;
                sA[threadIdx.x].q       = temp.w;
                sA[threadIdx.x].sig     = temp1.x;
                sA[threadIdx.x].eps     = temp1.y;
            }
            sA[threadIdx.x].fx      = 0.0f;
            sA[threadIdx.x].fy      = 0.0f;
            sA[threadIdx.x].fz      = 0.0f;
            apos.w                 *= cSim.epsfac;
            if (!bExclusionFlag)
            {
267
268
269
#ifdef USE_CUTOFF
                unsigned int flags = cSim.pInteractionFlag[pos];
                if (flags == 0)
270
                {
271
272
273
274
275
276
277
278
279
280
281
282
                    // No interactions in this block.
                }
                else if (flags == 0xFFFFFFFF)
#endif
                {
                    // Compute all interactions within this block.

                    for (unsigned int j = 0; j < GRID; j++)
                    {
                        dx              = psA[tj].x - apos.x;
                        dy              = psA[tj].y - apos.y;
                        dz              = psA[tj].z - apos.z;
283
#ifdef USE_PERIODIC
284
285
286
                        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;
287
#endif
288
289
290
291
292
293
294
295
                        r2              = dx * dx + dy * dy + dz * dz;
                        invR            = 1.0f / sqrt(r2);
                        sig             = a.x + psA[tj].sig;
                        sig2            = invR * sig;
                        sig2           *= sig2;
                        sig6            = sig2 * sig2 * sig2;
                        eps             = a.y * psA[tj].eps;
                        dEdR            = eps * (12.0f * sig6 - 6.0f) * sig6;
296
297
			/* E */
			CDLJ_energy     = eps * (sig6 - 1.0f) * sig6;
298
#ifdef USE_CUTOFF
299
    #ifdef USE_EWALD
300
301
302
303
304
                        float r         = sqrt(r2);
                        float alphaR    = cSim.alphaEwald * r;
                        dEdR           += apos.w * psA[tj].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
                        /* E */
                        CDLJ_energy    += apos.w * psA[tj].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
305
    #else
306
                        dEdR           += apos.w * psA[tj].q * (invR - 2.0f * cSim.reactionFieldK * r2);
307
308
			/* E */
                        CDLJ_energy    += apos.w * psA[tj].q * (invR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
309
    #endif
310
#else
311
                        dEdR           += apos.w * psA[tj].q * invR;
312
313
                        /* E */
                        CDLJ_energy    += apos.w * psA[tj].q * invR;
314
#endif
315
                        dEdR           *= invR * invR;
316
#ifdef USE_CUTOFF
317
318
319
                        if (r2 > cSim.nonbondedCutoffSqr)
                        {
                            dEdR = 0.0f;
320
321
			    /* E */
       			    CDLJ_energy = 0.0f;
322
323
                        }
#endif
324
325
			/* E */
			energy         += CDLJ_energy;
326
327
328
329
330
331
332
333
334
335
                        dx             *= dEdR;
                        dy             *= dEdR;
                        dz             *= dEdR;
                        af.x           -= dx;
                        af.y           -= dy;
                        af.z           -= dz;
                        psA[tj].fx     += dx;
                        psA[tj].fy     += dy;
                        psA[tj].fz     += dz;
                        tj              = (tj + 1) & (GRID - 1);
336
                    }
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
                }
#ifdef USE_CUTOFF
                else
                {
                    // Compute only a subset of the interactions in this block.

                    for (unsigned int j = 0; j < GRID; j++)
                    {
                        if ((flags&(1<<j)) != 0)
                        {
                            dx              = psA[j].x - apos.x;
                            dy              = psA[j].y - apos.y;
                            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;
354
#endif
355
356
357
358
359
360
361
362
                            r2              = dx * dx + dy * dy + dz * dz;
                            invR            = 1.0f / sqrt(r2);
                            sig             = a.x + psA[j].sig;
                            sig2            = invR * sig;
                            sig2           *= sig2;
                            sig6            = sig2 * sig2 * sig2;
                            eps             = a.y * psA[j].eps;
                            dEdR            = eps * (12.0f * sig6 - 6.0f) * sig6;
363
364
			    /* E */
			    CDLJ_energy     = eps * (sig6 - 1.0f) * sig6;
365
#ifdef USE_CUTOFF
366
    #ifdef USE_EWALD
367
368
369
370
                            float r         = sqrt(r2);
                            float alphaR    = cSim.alphaEwald * r;
                            dEdR           += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
                            CDLJ_energy    += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
371
    #else
372
                            dEdR           += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2);
373
374
                            /* E */
                            CDLJ_energy    += apos.w * psA[j].q * (invR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
375
    #endif
376
377
#else
                            dEdR           += apos.w * psA[j].q * invR;
378
379
                            /* E */
                            CDLJ_energy    += apos.w * psA[j].q * invR;
380
381
382
383
384
385
#endif
                            dEdR           *= invR * invR;
#ifdef USE_CUTOFF
                            if (r2 > cSim.nonbondedCutoffSqr)
                            {
                                dEdR = 0.0f;
386
387
				/* E */
				CDLJ_energy = 0.0f;
388
389
                            }
#endif
390
391
			    /* E */
			    energy         += CDLJ_energy;
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
                            dx             *= dEdR;
                            dy             *= dEdR;
                            dz             *= dEdR;
                            af.x           -= dx;
                            af.y           -= dy;
                            af.z           -= dz;
                            tempBuffer[threadIdx.x].x = dx;
                            tempBuffer[threadIdx.x].y = dy;
                            tempBuffer[threadIdx.x].z = dz;

                            // Sum the forces on atom j.

                            if (tgx % 2 == 0)
                            {
                                tempBuffer[threadIdx.x].x += tempBuffer[threadIdx.x+1].x;
                                tempBuffer[threadIdx.x].y += tempBuffer[threadIdx.x+1].y;
                                tempBuffer[threadIdx.x].z += tempBuffer[threadIdx.x+1].z;
                            }
                            if (tgx % 4 == 0)
                            {
                                tempBuffer[threadIdx.x].x += tempBuffer[threadIdx.x+2].x;
                                tempBuffer[threadIdx.x].y += tempBuffer[threadIdx.x+2].y;
                                tempBuffer[threadIdx.x].z += tempBuffer[threadIdx.x+2].z;
                            }
                            if (tgx % 8 == 0)
                            {
                                tempBuffer[threadIdx.x].x += tempBuffer[threadIdx.x+4].x;
                                tempBuffer[threadIdx.x].y += tempBuffer[threadIdx.x+4].y;
                                tempBuffer[threadIdx.x].z += tempBuffer[threadIdx.x+4].z;
                            }
                            if (tgx % 16 == 0)
                            {
                                tempBuffer[threadIdx.x].x += tempBuffer[threadIdx.x+8].x;
                                tempBuffer[threadIdx.x].y += tempBuffer[threadIdx.x+8].y;
                                tempBuffer[threadIdx.x].z += tempBuffer[threadIdx.x+8].z;
                            }
                            if (tgx == 0)
                            {
                                psA[j].fx += tempBuffer[threadIdx.x].x + tempBuffer[threadIdx.x+16].x;
                                psA[j].fy += tempBuffer[threadIdx.x].y + tempBuffer[threadIdx.x+16].y;
                                psA[j].fz += tempBuffer[threadIdx.x].z + tempBuffer[threadIdx.x+16].z;
                            }
                        }
                    }
436
                }
437
#endif
438
439
440
441
            }
            else  // bExclusion
            {
                // Read fixed atom data into registers and GRF
442
443
                unsigned int xi   = x>>GRIDBITS;
                unsigned int yi   = y>>GRIDBITS;
444
                unsigned int cell          = xi+yi*cSim.paddedNumberOfAtoms/GRID-yi*(yi+1)/2;
445
446
                unsigned int excl = cSim.pExclusion[cSim.pExclusionIndex[cell]+tgx];
                excl              = (excl >> tgx) | (excl << (GRID - tgx));
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
                for (unsigned int j = 0; j < GRID; j++)
                {
                    dx              = psA[tj].x - apos.x;
                    dy              = psA[tj].y - apos.y;
                    dz              = psA[tj].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
                    r2              = dx * dx + dy * dy + dz * dz;
                    invR            = 1.0f / sqrt(r2);
                    sig             = a.x + psA[tj].sig;
                    sig2            = invR * sig;
                    sig2           *= sig2;
                    sig6            = sig2 * sig2 * sig2;
                    eps             = a.y * psA[tj].eps;
                    dEdR            = eps * (12.0f * sig6 - 6.0f) * sig6;
465
466
		    /* E */
		    CDLJ_energy     = eps * (sig6 - 1.0f) * sig6;
467
#ifdef USE_CUTOFF
468
    #ifdef USE_EWALD
469
                    float r         = sqrt(r2);
470
                    float alphaR    = cSim.alphaEwald * r;
471
472
473
                    dEdR           += apos.w * psA[tj].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
                    /* E */
                    CDLJ_energy    += apos.w * psA[tj].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
474
    #else
475
                    dEdR           += apos.w * psA[tj].q * (invR - 2.0f * cSim.reactionFieldK * r2);
476
477
                    /* E */
		    CDLJ_energy    += apos.w * psA[tj].q * (invR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
478
    #endif
479
480
#else
                    dEdR           += apos.w * psA[tj].q * invR;
481
482
                    /* E */
                    CDLJ_energy    += apos.w * psA[tj].q * invR;
483
484
485
486
487
488
489
490
#endif
                    dEdR           *= invR * invR;
#ifdef USE_CUTOFF
                    if (!(excl & 0x1) || r2 > cSim.nonbondedCutoffSqr)
#else
                    if (!(excl & 0x1))
#endif
                    {
491
492
493
                        dEdR = 0.0f;			
                        /* E */
			CDLJ_energy  = 0.0f;
494
                    }
495
496
		    /* E */
		    energy         += CDLJ_energy;
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
                    dx             *= dEdR;
                    dy             *= dEdR;
                    dz             *= dEdR;
                    af.x           -= dx;
                    af.y           -= dy;
                    af.z           -= dz;
                    psA[tj].fx     += dx;
                    psA[tj].fy     += dy;
                    psA[tj].fz     += dz;
                    excl          >>= 1;
                    tj              = (tj + 1) & (GRID - 1);
                }
            }

            // Write results
            float4 of;
#ifdef USE_OUTPUT_BUFFER_PER_WARP
514
            unsigned int offset                          = x + tgx + warp*cSim.stride;
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
            of                                  = cSim.pForce4a[offset];
            of.x                               += af.x;
            of.y                               += af.y;
            of.z                               += af.z;
            cSim.pForce4a[offset]               = of;
            offset                              = y + tgx + warp*cSim.stride;
            of                                  = cSim.pForce4a[offset];
            of.x                               += sA[threadIdx.x].fx;
            of.y                               += sA[threadIdx.x].fy;
            of.z                               += sA[threadIdx.x].fz;
            cSim.pForce4a[offset]               = of;
#else
            of.x                                = af.x;
            of.y                                = af.y;
            of.z                                = af.z;
            of.w                                = 0.0f;
531
            unsigned int offset                          = x + tgx + (y >> GRIDBITS) * cSim.stride;
532
533
534
535
536
537
538
539
540
541
542
543
            cSim.pForce4a[offset]               = of;
            of.x                                = sA[threadIdx.x].fx;
            of.y                                = sA[threadIdx.x].fy;
            of.z                                = sA[threadIdx.x].fz;
            offset                              = y + tgx + (x >> GRIDBITS) * cSim.stride;
            cSim.pForce4a[offset]               = of;
#endif
            lasty = y;
        }

        pos++;
    }
544
    gpu_energies[energyIndex] = energy;
545
}