kCalculateAmoebaCudaPmeFixedEField.cu 23.4 KB
Newer Older
1
2
3
4
5
6
7
8

//-----------------------------------------------------------------------------------------

//-----------------------------------------------------------------------------------------

#include "amoebaCudaKernels.h"
#include "kCalculateAmoebaCudaUtilities.h"

Mark Friedrichs's avatar
Mark Friedrichs committed
9
//#define AMOEBA_DEBUG
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

static __constant__ cudaGmxSimulation cSim;
static __constant__ cudaAmoebaGmxSimulation cAmoebaSim;

void SetCalculateAmoebaCudaPmeFixedEFieldSim(amoebaGpuContext amoebaGpu)
{
    cudaError_t status;
    gpuContext gpu = amoebaGpu->gpuContext;
    status         = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));    
    RTERROR(status, "SetCalculateAmoebaCudaPmeFixedEFieldSim: cudaMemcpyToSymbol: SetSim copy to cSim failed");
    status         = cudaMemcpyToSymbol(cAmoebaSim, &amoebaGpu->amoebaSim, sizeof(cudaAmoebaGmxSimulation));    
    RTERROR(status, "SetCalculateAmoebaCudaPmeFixedEFieldSim: cudaMemcpyToSymbol: SetSim copy to cAmoebaSim failed");
}

void GetCalculateAmoebaCudaPmeFixedEFieldSim(amoebaGpuContext amoebaGpu)
{
    cudaError_t status;
    gpuContext gpu = amoebaGpu->gpuContext;
    status = cudaMemcpyFromSymbol(&gpu->sim, cSim, sizeof(cudaGmxSimulation));    
    RTERROR(status, "GetCalculateAmoebaCudaPmeFixedEFieldSim: cudaMemcpyFromSymbol: SetSim copy from cSim failed");
    status = cudaMemcpyFromSymbol(&amoebaGpu->amoebaSim, cAmoebaSim, sizeof(cudaAmoebaGmxSimulation));         
    RTERROR(status, "GetCalculateAmoebaCudaPmeFixedEFieldSim: cudaMemcpyFromSymbol: SetSim copy from cAmoebaSim failed");
}

__global__ 
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
#elif (__CUDA_ARCH__ >= 130)
__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
#endif
Mark Friedrichs's avatar
Mark Friedrichs committed
42
static void kReducePmeEFieldPolar_kernel( unsigned int fieldComponents, unsigned int outputBuffers, float* EFieldReciprocal,  float* fieldIn, float* fieldOut )
43
44
45
46
47
{
    unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;

    // Reduce field

48
    const float term = (4.0f/3.0f)*(cSim.alphaEwald*cSim.alphaEwald*cSim.alphaEwald)/cAmoebaSim.sqrtPi;
Mark Friedrichs's avatar
Mark Friedrichs committed
49
    //const float term = 0.0f;
50
51
52
53
54
    while (pos < fieldComponents)
    {   

        // self-term included here

Mark Friedrichs's avatar
Mark Friedrichs committed
55
        float totalField = EFieldReciprocal[pos] + term*cAmoebaSim.pLabFrameDipole[pos];
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76

        float* pFt       = fieldIn + pos;
        unsigned int i   = outputBuffers;
        while (i >= 4)
        {   
            totalField += pFt[0] + pFt[fieldComponents] + pFt[2*fieldComponents] + pFt[3*fieldComponents];
            pFt        += fieldComponents*4;
            i          -= 4;
        }   

        if (i >= 2)
        {   
            totalField += pFt[0] + pFt[fieldComponents];
            pFt        += fieldComponents*2;
            i          -= 2;
        }   

        if (i > 0)
        {   
            totalField += pFt[0];
        }   
Mark Friedrichs's avatar
Mark Friedrichs committed
77

78
79
80
81
82
        fieldOut[pos]   = totalField;
        pos            += gridDim.x * blockDim.x;
    }   
}

Mark Friedrichs's avatar
Mark Friedrichs committed
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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
__global__ 
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
#elif (__CUDA_ARCH__ >= 130)
__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
#endif
static void kReducePmeEField_kernel( unsigned int fieldComponents, unsigned int outputBuffers,  float* fieldIn, float* fieldOut )
{
    unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;

    // Reduce field

    const float term = (4.0f/3.0f)*(cSim.alphaEwald*cSim.alphaEwald*cSim.alphaEwald)/cAmoebaSim.sqrtPi;
    while (pos < fieldComponents)
    {   

        // self-term included here

        float totalField = term*cAmoebaSim.pLabFrameDipole[pos];

        float* pFt       = fieldIn + pos;
        unsigned int i   = outputBuffers;
        while (i >= 4)
        {   
            totalField += pFt[0] + pFt[fieldComponents] + pFt[2*fieldComponents] + pFt[3*fieldComponents];
            pFt        += fieldComponents*4;
            i          -= 4;
        }   

        if (i >= 2)
        {   
            totalField += pFt[0] + pFt[fieldComponents];
            pFt        += fieldComponents*2;
            i          -= 2;
        }   

        if (i > 0)
        {   
            totalField += pFt[0];
        }   

        fieldOut[pos]  += totalField;
        pos            += gridDim.x * blockDim.x;
    }   
}

131
132
133
134
135
// reduce psWorkArray_3_1 -> EField
// reduce psWorkArray_3_2 -> EFieldPolar

static void kReducePmeDirectE_Fields(amoebaGpuContext amoebaGpu )
{
Mark Friedrichs's avatar
Mark Friedrichs committed
136
137
138
139
140
141

    // E_FieldPolar = E_Field (reciprocal) + E_FieldPolar (direct) + self

    kReducePmeEFieldPolar_kernel<<<amoebaGpu->nonbondBlocks, amoebaGpu->fieldReduceThreadsPerBlock>>>(
                                   amoebaGpu->paddedNumberOfAtoms*3, amoebaGpu->outputBuffers,
                                   amoebaGpu->psE_Field->_pDevStream[0], amoebaGpu->psWorkArray_3_2->_pDevStream[0], amoebaGpu->psE_FieldPolar->_pDevStream[0] );
142
143
    LAUNCHERROR("kReducePmeE_Fields1");

Mark Friedrichs's avatar
Mark Friedrichs committed
144
145
146
147
148
    // E_Field = E_Field (reciprocal) + E_Field (direct) + self

    kReducePmeEField_kernel<<<amoebaGpu->nonbondBlocks, amoebaGpu->fieldReduceThreadsPerBlock>>>(
                              amoebaGpu->paddedNumberOfAtoms*3, amoebaGpu->outputBuffers,
                              amoebaGpu->psWorkArray_3_1->_pDevStream[0], amoebaGpu->psE_Field->_pDevStream[0] );
149
150
151
152
153
154
    LAUNCHERROR("kReducePmeE_Fields2");
}

// file includes FixedFieldParticle struct definition/load/unload struct and body kernel for fixed E-field

#undef GK
Mark Friedrichs's avatar
Mark Friedrichs committed
155
156
#undef INCLUDE_FIXED_FIELD_BUFFERS
#define INCLUDE_FIXED_FIELD_BUFFERS
157
#include "kCalculateAmoebaCudaFixedFieldParticle.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
158
159
160
161
162
163
164
165
166
167
168
#undef INCLUDE_FIXED_FIELD_BUFFERS
__device__ void sumTempBuffer( FixedFieldParticle& atomI, FixedFieldParticle& atomJ ){
    atomI.tempBuffer[0]  += atomJ.tempBuffer[0];
    atomI.tempBuffer[1]  += atomJ.tempBuffer[1];
    atomI.tempBuffer[2]  += atomJ.tempBuffer[2];

    atomI.tempBufferP[0] += atomJ.tempBufferP[0];
    atomI.tempBufferP[1] += atomJ.tempBufferP[1];
    atomI.tempBufferP[2] += atomJ.tempBufferP[2];
}

169
__device__ void calculateFixedFieldRealSpacePairIxn_kernel( FixedFieldParticle& atomI, FixedFieldParticle& atomJ,
170
                                                            float dscale, float pscale, float4 fields[3]
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
#ifdef AMOEBA_DEBUG
                                                            , float4* pullBack
#endif

 ){

    // compute the real space portion of the Ewald summation
  
    float xr          = atomJ.x - atomI.x;
    float yr          = atomJ.y - atomI.y;
    float zr          = atomJ.z - atomI.z;

    // periodic boundary conditions

    xr               -= floor(xr*cSim.invPeriodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
    yr               -= floor(yr*cSim.invPeriodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
    zr               -= floor(zr*cSim.invPeriodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;

Mark Friedrichs's avatar
Mark Friedrichs committed
189
    float r2          = xr*xr + yr*yr + zr*zr;
190
191
192
193
    float r           = sqrtf(r2);

    // calculate the error function damping terms

194
    float ralpha      = cSim.alphaEwald*r;
195

196
    float bn0             = erfc(ralpha)/r;
197
198
    float alsq2       = 2.0f*cSim.alphaEwald*cSim.alphaEwald;
    float alsq2n      = 1.0f/(cAmoebaSim.sqrtPi*cSim.alphaEwald);
199
200
    float exp2a       = exp(-(ralpha*ralpha));
    alsq2n           *= alsq2;
201
    float bn1             = (bn0+alsq2n*exp2a)/r2;
202
203

    alsq2n           *= alsq2;
204
    float bn2             = (3.0f*bn1+alsq2n*exp2a)/r2;
205
206

    alsq2n           *= alsq2;
207
    float bn3             = (5.0f*bn2+alsq2n*exp2a)/r2;
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263

    // compute the error function scaled and unscaled terms

    float scale3      = 1.0f;
    float scale5      = 1.0f;
    float scale7      = 1.0f;
    float damp        = atomI.damp*atomJ.damp;
    if( damp != 0.0f ){

        float ratio  = (r/damp);
              ratio  = ratio*ratio*ratio;

        float pgamma = atomI.thole < atomJ.thole ? atomI.thole : atomJ.thole;

              damp   = -pgamma*ratio;

        if( damp > -50.0f) {
            float expdamp = exp(damp);
            scale3        = 1.0f - expdamp;
            scale5        = 1.0f - expdamp*(1.0f-damp);
            scale7        = 1.0f - expdamp*(1.0f-damp+(0.6f*damp*damp));
        }
    }
    float dsc3        = dscale*scale3;
    float dsc5        = dscale*scale5;
    float dsc7        = dscale*scale7;

    float psc3        = pscale*scale3;
    float psc5        = pscale*scale5;
    float psc7        = pscale*scale7;

    float r3          = (r*r2);
    float r5          = (r3*r2);
    float r7          = (r5*r2);
    float drr3        = (1.0f-dsc3)/r3;
    float drr5        = 3.0f * (1.0f-dsc5)/r5;
    float drr7        = 15.0f * (1.0f-dsc7)/r7;

    float prr3        = (1.0f-psc3) / r3;
    float prr5        = 3.0f *(1.0f-psc5)/r5;
    float prr7        = 15.0f*(1.0f-psc7)/r7;

    float dir         = atomI.labFrameDipole_X*xr + atomI.labFrameDipole_Y*yr + atomI.labFrameDipole_Z*zr;

    float qix         = atomI.labFrameQuadrupole_XX*xr + atomI.labFrameQuadrupole_XY*yr + atomI.labFrameQuadrupole_XZ*zr;
    float qiy         = atomI.labFrameQuadrupole_XY*xr + atomI.labFrameQuadrupole_YY*yr + atomI.labFrameQuadrupole_YZ*zr;
    float qiz         = atomI.labFrameQuadrupole_XZ*xr + atomI.labFrameQuadrupole_YZ*yr + atomI.labFrameQuadrupole_ZZ*zr;

    float qir         = qix*xr + qiy*yr + qiz*zr;

    float dkr         = atomJ.labFrameDipole_X*xr + atomJ.labFrameDipole_Y*yr + atomJ.labFrameDipole_Z*zr;
    float qkx         = atomJ.labFrameQuadrupole_XX*xr + atomJ.labFrameQuadrupole_XY*yr + atomJ.labFrameQuadrupole_XZ*zr;
    float qky         = atomJ.labFrameQuadrupole_XY*xr + atomJ.labFrameQuadrupole_YY*yr + atomJ.labFrameQuadrupole_YZ*zr;
    float qkz         = atomJ.labFrameQuadrupole_XZ*xr + atomJ.labFrameQuadrupole_YZ*yr + atomJ.labFrameQuadrupole_ZZ*zr;
    float qkr         = qkx*xr + qky*yr + qkz*zr;

264
265
    float fim0            = -xr*(bn1*atomJ.q-bn2*dkr+bn3*qkr)
                         - bn1*atomJ.labFrameDipole_X + 2.0f*bn2*qkx;
266

267
268
    float fim1            = -yr*(bn1*atomJ.q-bn2*dkr+bn3*qkr)
                         - bn1*atomJ.labFrameDipole_Y + 2.0f*bn2*qky;
269

270
271
    float fim2            = -zr*(bn1*atomJ.q-bn2*dkr+bn3*qkr)
                         - bn1*atomJ.labFrameDipole_Z + 2.0f*bn2*qkz;
272

273
274
    float fkm0            = xr*(bn1*atomI.q+bn2*dir+bn3*qir)
                         - bn1*atomI.labFrameDipole_X - 2.0f*bn2*qix;
275

276
277
    float fkm1            = yr*(bn1*atomI.q+bn2*dir+bn3*qir)
                         - bn1*atomI.labFrameDipole_Y - 2.0f*bn2*qiy;
278

279
280
    float fkm2            = zr*(bn1*atomI.q+bn2*dir+bn3*qir)
                         - bn1*atomI.labFrameDipole_Z - 2.0f*bn2*qiz;
281

282
    float fid0            = -xr*(drr3*atomJ.q-drr5*dkr+drr7*qkr)
283
284
                         - drr3*atomJ.labFrameDipole_X + 2.0f*drr5*qkx;

285
    float fid1            = -yr*(drr3*atomJ.q-drr5*dkr+drr7*qkr)
286
287
                         - drr3*atomJ.labFrameDipole_Y + 2.0f*drr5*qky;

288
    float fid2            = -zr*(drr3*atomJ.q-drr5*dkr+drr7*qkr)
289
290
                         - drr3*atomJ.labFrameDipole_Z + 2.0f*drr5*qkz;

291
    float fkd0            = xr*(drr3*atomI.q+drr5*dir+drr7*qir)
292
293
                         - drr3*atomI.labFrameDipole_X - 2.0f*drr5*qix;

294
    float fkd1            = yr*(drr3*atomI.q+drr5*dir+drr7*qir)
295
296
                         - drr3*atomI.labFrameDipole_Y - 2.0f*drr5*qiy;

297
    float fkd2            = zr*(drr3*atomI.q+drr5*dir+drr7*qir)
298
299
                         - drr3*atomI.labFrameDipole_Z - 2.0f*drr5*qiz;

300
    float fip0            = -xr*(prr3*atomJ.q-prr5*dkr+prr7*qkr)
301
302
                         - prr3*atomJ.labFrameDipole_X + 2.0f*prr5*qkx;

303
    float fip1            = -yr*(prr3*atomJ.q-prr5*dkr+prr7*qkr)
304
305
                         - prr3*atomJ.labFrameDipole_Y + 2.0f*prr5*qky;

306
    float fip2            = -zr*(prr3*atomJ.q-prr5*dkr+prr7*qkr)
307
308
                         - prr3*atomJ.labFrameDipole_Z + 2.0f*prr5*qkz;

309
    float fkp0            = xr*(prr3*atomI.q+prr5*dir+prr7*qir)
310
311
                         - prr3*atomI.labFrameDipole_X - 2.0f*prr5*qix;

312
    float fkp1            = yr*(prr3*atomI.q+prr5*dir+prr7*qir)
313
314
                         - prr3*atomI.labFrameDipole_Y - 2.0f*prr5*qiy;

315
    float fkp2            = zr*(prr3*atomI.q+prr5*dir+prr7*qir)
316
317
318
319
                         - prr3*atomI.labFrameDipole_Z - 2.0f*prr5*qiz;
  
    // increment the field at each site due to this interaction

Mark Friedrichs's avatar
Mark Friedrichs committed
320
    if( r2 <= cSim.nonbondedCutoffSqr ){
321

322
323
324
        fields[0].x       = fim0 - fid0;
        fields[1].x       = fim1 - fid1;
        fields[2].x       = fim2 - fid2;
Mark Friedrichs's avatar
Mark Friedrichs committed
325

326
327
328
        fields[0].y       = fkm0 - fkd0;
        fields[1].y       = fkm1 - fkd1;
        fields[2].y       = fkm2 - fkd2;
Mark Friedrichs's avatar
Mark Friedrichs committed
329

330
331
332
        fields[0].z       = fim0 - fip0;
        fields[1].z       = fim1 - fip1;
        fields[2].z       = fim2 - fip2;
Mark Friedrichs's avatar
Mark Friedrichs committed
333

334
335
336
        fields[0].w       = fkm0 - fkp0;
        fields[1].w       = fkm1 - fkp1;
        fields[2].w       = fkm2 - fkp2;
337
338
339
 
    } else {

340
341
342
343
        fields[0].x       = 0.0f;
        fields[0].y       = 0.0f;
        fields[0].z       = 0.0f;
        fields[0].w       = 0.0f;
344
    
345
346
347
348
        fields[1].x       = 0.0f;
        fields[1].y       = 0.0f;
        fields[1].z       = 0.0f;
        fields[1].w       = 0.0f;
349
    
350
351
352
353
        fields[2].x       = 0.0f;
        fields[2].y       = 0.0f;
        fields[2].z       = 0.0f;
        fields[2].w       = 0.0f;
354
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
355

356
357
358
359
360
361
362
363
364
365
#ifdef AMOEBA_DEBUG
    pullBack[0].x = xr;
    pullBack[0].y = yr;
    pullBack[0].z = zr;
    pullBack[0].w = r2;

    pullBack[1].x = atomJ.x - atomI.x;
    pullBack[1].y = atomJ.y - atomI.y;
    pullBack[1].z = atomJ.z - atomI.z;
    pullBack[1].w = (atomJ.x - atomI.x)*(atomJ.x - atomI.x) + (atomJ.y - atomI.y)*(atomJ.y - atomI.y)+ (atomJ.z - atomI.z)*(atomJ.z - atomI.z);
Mark Friedrichs's avatar
Mark Friedrichs committed
366
367
368
369
370
371

    pullBack[2].x = scale3;
    pullBack[2].y = scale5;
    pullBack[2].z = scale7;
    pullBack[2].w = -1.0;

372
373
374
375
376
377
378
379
380
381
382
383
#endif
}

// Include versions of the kernels for N^2 calculations.

#define METHOD_NAME(a, b) a##N2##b
#include "kCalculateAmoebaCudaPmeFixedEField.h"
#define USE_OUTPUT_BUFFER_PER_WARP
#undef METHOD_NAME
#define METHOD_NAME(a, b) a##N2ByWarp##b
#include "kCalculateAmoebaCudaPmeFixedEField.h"

Mark Friedrichs's avatar
Mark Friedrichs committed
384
385
386
387
388
389
390
391
392
/**---------------------------------------------------------------------------------------

   Report whether a number is a nan or infinity

   @param number               number to test
   @return 1 if number is  nan or infinity; else return 0

   --------------------------------------------------------------------------------------- */

Mark Friedrichs's avatar
Mark Friedrichs committed
393
#ifdef AMOEBA_DEBUG
Mark Friedrichs's avatar
Mark Friedrichs committed
394
395
396
static int isNanOrInfinity( double number ){
    return (number != number || number == std::numeric_limits<double>::infinity() || number == -std::numeric_limits<double>::infinity()) ? 1 : 0; 
}
Mark Friedrichs's avatar
Mark Friedrichs committed
397
#endif
Mark Friedrichs's avatar
Mark Friedrichs committed
398

399
400
401
402
403
404
405
406
/**---------------------------------------------------------------------------------------

   Compute fixed electric field using PME

   @param amoebaGpu        amoebaGpu context

   --------------------------------------------------------------------------------------- */

Mark Friedrichs's avatar
Mark Friedrichs committed
407
static void cudaComputeAmoebaPmeDirectFixedEField( amoebaGpuContext amoebaGpu )
408
409
{
  
Mark Friedrichs's avatar
Mark Friedrichs committed
410
411
    static unsigned int threadsPerBlock  = 0;
    gpuContext gpu                       = amoebaGpu->gpuContext;
412
413
414
415
416
417
418
419
420
421
422

#ifdef AMOEBA_DEBUG
    static const char* methodName = "computeCudaAmoebaPmeFixedEField";
    if( amoebaGpu->log ){
        (void) fprintf( amoebaGpu->log, "\n%s\n", methodName ); (void) fflush( amoebaGpu->log );
    }
    int paddedNumberOfAtoms                    = amoebaGpu->gpuContext->sim.paddedNumberOfAtoms;

    // N2 debug array

    CUDAStream<float4>* debugArray             = new CUDAStream<float4>(paddedNumberOfAtoms*paddedNumberOfAtoms, 1, "DebugArray");
Mark Friedrichs's avatar
Mark Friedrichs committed
423
    memset( debugArray->_pSysStream[0], 0, sizeof( float )*4*paddedNumberOfAtoms*paddedNumberOfAtoms);
424
425
426
427
    debugArray->Upload();

    // print intermediate results for the targetAtom 

Mark Friedrichs's avatar
Mark Friedrichs committed
428
    unsigned int targetAtom  = 1280;
Mark Friedrichs's avatar
Mark Friedrichs committed
429
#endif
Mark Friedrichs's avatar
Mark Friedrichs committed
430

431
432
    kClearFields_3( amoebaGpu, 2 );

Mark Friedrichs's avatar
Mark Friedrichs committed
433
434
435
436
437
438
439
    // on first pass, set threads/block

    if( threadsPerBlock == 0 ){ 
        unsigned int maxThreads;
        if (gpu->sm_version >= SM_20)
            maxThreads = 384; 
        else if (gpu->sm_version >= SM_12)
440
            maxThreads = 192;
Mark Friedrichs's avatar
Mark Friedrichs committed
441
442
443
444
445
        else
            maxThreads = 64;
        threadsPerBlock = std::min(getThreadsPerBlock(amoebaGpu, sizeof(FixedFieldParticle)), maxThreads);
    }    

446
    if (gpu->bOutputBufferPerWarp){
Mark Friedrichs's avatar
Mark Friedrichs committed
447
448
        kCalculateAmoebaPmeDirectFixedE_FieldN2ByWarp_kernel<<<amoebaGpu->nonbondBlocks, threadsPerBlock, sizeof(FixedFieldParticle)*threadsPerBlock>>>(
                                                                           gpu->sim.pInteractingWorkUnit,
449
450
451
452
453
454
455
456
                                                                           amoebaGpu->psWorkArray_3_1->_pDevStream[0],
#ifdef AMOEBA_DEBUG
                                                                           amoebaGpu->psWorkArray_3_2->_pDevStream[0],
                                                                           debugArray->_pDevStream[0], targetAtom );
#else
                                                                           amoebaGpu->psWorkArray_3_2->_pDevStream[0] );
#endif
    } else {
Mark Friedrichs's avatar
Mark Friedrichs committed
457
458
        kCalculateAmoebaPmeDirectFixedE_FieldN2_kernel<<<amoebaGpu->nonbondBlocks, threadsPerBlock, sizeof(FixedFieldParticle)*threadsPerBlock>>>(
                                                                           gpu->sim.pInteractingWorkUnit,
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
                                                                           amoebaGpu->psWorkArray_3_1->_pDevStream[0],
#ifdef AMOEBA_DEBUG
                                                                           amoebaGpu->psWorkArray_3_2->_pDevStream[0],
                                                                           debugArray->_pDevStream[0], targetAtom );
#else
                                                                           amoebaGpu->psWorkArray_3_2->_pDevStream[0] );
#endif
    }
    LAUNCHERROR("kCalculateAmoebaPmeDirectFixedE_Field_kernel");

    kReducePmeDirectE_Fields( amoebaGpu );

#ifdef AMOEBA_DEBUG
    if( amoebaGpu->log ){
        gpu->psInteractionCount->Download();
Mark Friedrichs's avatar
Mark Friedrichs committed
474
475
476
        (void) fprintf( amoebaGpu->log, "cudaComputeAmoebaPmeDirectFixedEField:  threadsPerBlock=%u getThreadsPerBlock=%d sizeof=%u shrd=%u\n", 
                        threadsPerBlock, getThreadsPerBlock(amoebaGpu, sizeof(FixedFieldParticle)+sizeof(float3)),
                        (sizeof(FixedFieldParticle)+sizeof(float3)), (sizeof(FixedFieldParticle)+sizeof(float3))*threadsPerBlock );
Mark Friedrichs's avatar
Mark Friedrichs committed
477
        (void) fprintf( amoebaGpu->log, "AmoebaN2Forces_kernel numBlocks=%u numThreads=%u bufferPerWarp=%u atm=%u shrd=%u Ebuf=%u ixnCt=%u workUnits=%u warp=%d\n",
Mark Friedrichs's avatar
Mark Friedrichs committed
478
                        amoebaGpu->nonbondBlocks, threadsPerBlock, amoebaGpu->bOutputBufferPerWarp,
Mark Friedrichs's avatar
Mark Friedrichs committed
479
480
                        sizeof(FixedFieldParticle), sizeof(FixedFieldParticle)*amoebaGpu->nonbondThreadsPerBlock, amoebaGpu->energyOutputBuffers, 
                        (*gpu->psInteractionCount)[0], gpu->sim.workUnits, gpu->bOutputBufferPerWarp );
481
        (void) fflush( amoebaGpu->log );
Mark Friedrichs's avatar
Mark Friedrichs committed
482
483
/*
        (void) fprintf( amoebaGpu->log, "Out WorkArray_3_[1,2]  paddedNumberOfAtoms=%d\n",  amoebaGpu->paddedNumberOfAtoms, amoebaGpu->outputBuffers );
484
485
        amoebaGpu->psWorkArray_3_1->Download();
        amoebaGpu->psWorkArray_3_2->Download();
Mark Friedrichs's avatar
Mark Friedrichs committed
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
        for( int ii = 0; ii < amoebaGpu->paddedNumberOfAtoms; ii++ ){
           (void) fprintf( amoebaGpu->log, "%5d ", ii); 

            int indexOffset     = ii*3;

           // buffer 1

           (void) fprintf( amoebaGpu->log,"WArry1[%16.9e %16.9e %16.9e] ",
                           amoebaGpu->psWorkArray_3_1->_pSysStream[0][indexOffset],
                           amoebaGpu->psWorkArray_3_1->_pSysStream[0][indexOffset+1],
                           amoebaGpu->psWorkArray_3_1->_pSysStream[0][indexOffset+2] );
   
           // buffer 2

           (void) fprintf( amoebaGpu->log,"WArry2[%16.9e %16.9e %16.9e] ",
                           amoebaGpu->psWorkArray_3_2->_pSysStream[0][indexOffset],
                           amoebaGpu->psWorkArray_3_2->_pSysStream[0][indexOffset+1],
                           amoebaGpu->psWorkArray_3_2->_pSysStream[0][indexOffset+2] );

           (void) fprintf( amoebaGpu->log,"\n" );
           if( ii == maxPrint && (gpu->natoms - maxPrint) > ii ){
                ii = gpu->natoms - maxPrint;
           }
        }
        (void) fflush( amoebaGpu->log );
Mark Friedrichs's avatar
Mark Friedrichs committed
511
512
513
*/
        amoebaGpu->psE_Field->Download();
        amoebaGpu->psE_FieldPolar->Download();
Mark Friedrichs's avatar
Mark Friedrichs committed
514
515
        (void) fprintf( amoebaGpu->log,"E-field (includes self term)" );
        int maxPrint             = 3002;
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
        for( int ii = 0; ii < gpu->natoms; ii++ ){
           (void) fprintf( amoebaGpu->log, "%5d ", ii); 

            int indexOffset     = ii*3;

           // E_Field

           (void) fprintf( amoebaGpu->log,"E[%16.9e %16.9e %16.9e] ",
                           amoebaGpu->psE_Field->_pSysStream[0][indexOffset],
                           amoebaGpu->psE_Field->_pSysStream[0][indexOffset+1],
                           amoebaGpu->psE_Field->_pSysStream[0][indexOffset+2] );
   
           // E_Field polar

           (void) fprintf( amoebaGpu->log,"Epol[%16.9e %16.9e %16.9e] ",
                           amoebaGpu->psE_FieldPolar->_pSysStream[0][indexOffset],
                           amoebaGpu->psE_FieldPolar->_pSysStream[0][indexOffset+1],
                           amoebaGpu->psE_FieldPolar->_pSysStream[0][indexOffset+2] );

           (void) fprintf( amoebaGpu->log,"\n" );
           if( ii == maxPrint && (gpu->natoms - maxPrint) > ii ){
                ii = gpu->natoms - maxPrint;
           }
        }
        (void) fflush( amoebaGpu->log );
        (void) fprintf( amoebaGpu->log, "EFields End\n" );

        (void) fprintf( amoebaGpu->log, "DebugQ\n" );
        debugArray->Download();

        int paddedNumberOfAtoms                    = amoebaGpu->gpuContext->sim.paddedNumberOfAtoms;
Mark Friedrichs's avatar
Mark Friedrichs committed
547
        amoebaGpu->gpuContext->psPosq4->Download();
548
549
        for( int jj = 0; jj < gpu->natoms; jj++ ){
            int debugIndex = jj;
Mark Friedrichs's avatar
Mark Friedrichs committed
550
if( fabs(debugArray->_pSysStream[0][jj+3*paddedNumberOfAtoms].x) > 0.0 ){
551
            (void) fprintf( amoebaGpu->log,"%5d PmeFixedEField\n", jj );
Mark Friedrichs's avatar
Mark Friedrichs committed
552
            for( int kk = 0; kk < 7; kk++ ){
553
554
555
556
557
558
                (void) fprintf( amoebaGpu->log,"[%16.9e %16.9e %16.9e %16.9e]\n",
                                debugArray->_pSysStream[0][debugIndex].x, debugArray->_pSysStream[0][debugIndex].y,
                                debugArray->_pSysStream[0][debugIndex].z, debugArray->_pSysStream[0][debugIndex].w );
                debugIndex += paddedNumberOfAtoms;
            }
            (void) fprintf( amoebaGpu->log,"\n" );
Mark Friedrichs's avatar
Mark Friedrichs committed
559
}
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577

        }

        // write results to file

        if( 1 ){
            std::vector<int> fileId;
            //fileId.push_back( 0 );
            VectorOfDoubleVectors outputVector;
            cudaLoadCudaFloat4Array( gpu->natoms, 3, gpu->psPosq4,              outputVector );
            cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psE_Field,      outputVector );
            cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psE_FieldPolar, outputVector);
            cudaWriteVectorOfDoubleVectorsToFile( "CudaEField", fileId, outputVector );
         }
         delete debugArray;
    }
#endif

Mark Friedrichs's avatar
Mark Friedrichs committed
578
        if( 0 ){
579
580
581
582
583
584
585
586
587
588
            std::vector<int> fileId;
            fileId.push_back( 0 );
            VectorOfDoubleVectors outputVector;
            cudaLoadCudaFloat4Array( gpu->natoms, 3, gpu->psPosq4,              outputVector );
            cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psE_Field,      outputVector );
            cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psE_FieldPolar, outputVector);
            cudaWriteVectorOfDoubleVectorsToFile( "CudaEField", fileId, outputVector );

         }
}
Mark Friedrichs's avatar
Mark Friedrichs committed
589
590
591

void cudaComputeAmoebaPmeFixedEField( amoebaGpuContext amoebaGpu )
{
592
    kCalculateAmoebaPMEFixedMultipoles( amoebaGpu );
Mark Friedrichs's avatar
Mark Friedrichs committed
593
    cudaComputeAmoebaPmeDirectFixedEField( amoebaGpu );
Mark Friedrichs's avatar
Mark Friedrichs committed
594
}