kCalculateAmoebaCudaPmeFixedEField.cu 23.8 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;
Peter Eastman's avatar
Peter Eastman committed
190
191
    if( r2 <= cSim.nonbondedCutoffSqr ){
        float r           = sqrtf(r2);
192

Peter Eastman's avatar
Peter Eastman committed
193
        // calculate the error function damping terms
194

Peter Eastman's avatar
Peter Eastman committed
195
        float ralpha      = cSim.alphaEwald*r;
196

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

Peter Eastman's avatar
Peter Eastman committed
204
205
        alsq2n           *= alsq2;
        float bn2             = (3.0f*bn1+alsq2n*exp2a)/r2;
206

Peter Eastman's avatar
Peter Eastman committed
207
208
        alsq2n           *= alsq2;
        float bn3             = (5.0f*bn2+alsq2n*exp2a)/r2;
209

Peter Eastman's avatar
Peter Eastman committed
210
        // compute the error function scaled and unscaled terms
211

Peter Eastman's avatar
Peter Eastman committed
212
213
214
215
216
        float scale3      = 1.0f;
        float scale5      = 1.0f;
        float scale7      = 1.0f;
        float damp        = atomI.damp*atomJ.damp;
        if( damp != 0.0f ){
217

Peter Eastman's avatar
Peter Eastman committed
218
219
            float ratio  = (r/damp);
                  ratio  = ratio*ratio*ratio;
220

Peter Eastman's avatar
Peter Eastman committed
221
            float pgamma = atomI.thole < atomJ.thole ? atomI.thole : atomJ.thole;
222

Peter Eastman's avatar
Peter Eastman committed
223
                  damp   = -pgamma*ratio;
224

Peter Eastman's avatar
Peter Eastman committed
225
226
227
228
229
230
            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));
            }
231
        }
Peter Eastman's avatar
Peter Eastman committed
232
233
234
        float dsc3        = dscale*scale3;
        float dsc5        = dscale*scale5;
        float dsc7        = dscale*scale7;
235

Peter Eastman's avatar
Peter Eastman committed
236
237
238
        float psc3        = pscale*scale3;
        float psc5        = pscale*scale5;
        float psc7        = pscale*scale7;
239

Peter Eastman's avatar
Peter Eastman committed
240
241
242
243
244
245
        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;
246

Peter Eastman's avatar
Peter Eastman committed
247
248
249
        float prr3        = (1.0f-psc3) / r3;
        float prr5        = 3.0f *(1.0f-psc5)/r5;
        float prr7        = 15.0f*(1.0f-psc7)/r7;
250

Peter Eastman's avatar
Peter Eastman committed
251
        float dir         = atomI.labFrameDipole_X*xr + atomI.labFrameDipole_Y*yr + atomI.labFrameDipole_Z*zr;
252

Peter Eastman's avatar
Peter Eastman committed
253
254
255
        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;
256

Peter Eastman's avatar
Peter Eastman committed
257
        float qir         = qix*xr + qiy*yr + qiz*zr;
258

Peter Eastman's avatar
Peter Eastman committed
259
260
261
262
263
        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

Peter Eastman's avatar
Peter Eastman committed
265
266
        float fim0            = -xr*(bn1*atomJ.q-bn2*dkr+bn3*qkr)
                             - bn1*atomJ.labFrameDipole_X + 2.0f*bn2*qkx;
267

Peter Eastman's avatar
Peter Eastman committed
268
269
        float fim1            = -yr*(bn1*atomJ.q-bn2*dkr+bn3*qkr)
                             - bn1*atomJ.labFrameDipole_Y + 2.0f*bn2*qky;
270

Peter Eastman's avatar
Peter Eastman committed
271
272
        float fim2            = -zr*(bn1*atomJ.q-bn2*dkr+bn3*qkr)
                             - bn1*atomJ.labFrameDipole_Z + 2.0f*bn2*qkz;
273

Peter Eastman's avatar
Peter Eastman committed
274
275
        float fkm0            = xr*(bn1*atomI.q+bn2*dir+bn3*qir)
                             - bn1*atomI.labFrameDipole_X - 2.0f*bn2*qix;
276

Peter Eastman's avatar
Peter Eastman committed
277
278
        float fkm1            = yr*(bn1*atomI.q+bn2*dir+bn3*qir)
                             - bn1*atomI.labFrameDipole_Y - 2.0f*bn2*qiy;
279

Peter Eastman's avatar
Peter Eastman committed
280
281
        float fkm2            = zr*(bn1*atomI.q+bn2*dir+bn3*qir)
                             - bn1*atomI.labFrameDipole_Z - 2.0f*bn2*qiz;
282

Peter Eastman's avatar
Peter Eastman committed
283
284
        float fid0            = -xr*(drr3*atomJ.q-drr5*dkr+drr7*qkr)
                             - drr3*atomJ.labFrameDipole_X + 2.0f*drr5*qkx;
285

Peter Eastman's avatar
Peter Eastman committed
286
287
        float fid1            = -yr*(drr3*atomJ.q-drr5*dkr+drr7*qkr)
                             - drr3*atomJ.labFrameDipole_Y + 2.0f*drr5*qky;
288

Peter Eastman's avatar
Peter Eastman committed
289
290
        float fid2            = -zr*(drr3*atomJ.q-drr5*dkr+drr7*qkr)
                             - drr3*atomJ.labFrameDipole_Z + 2.0f*drr5*qkz;
291

Peter Eastman's avatar
Peter Eastman committed
292
293
        float fkd0            = xr*(drr3*atomI.q+drr5*dir+drr7*qir)
                             - drr3*atomI.labFrameDipole_X - 2.0f*drr5*qix;
294

Peter Eastman's avatar
Peter Eastman committed
295
296
        float fkd1            = yr*(drr3*atomI.q+drr5*dir+drr7*qir)
                             - drr3*atomI.labFrameDipole_Y - 2.0f*drr5*qiy;
297

Peter Eastman's avatar
Peter Eastman committed
298
299
        float fkd2            = zr*(drr3*atomI.q+drr5*dir+drr7*qir)
                             - drr3*atomI.labFrameDipole_Z - 2.0f*drr5*qiz;
300

Peter Eastman's avatar
Peter Eastman committed
301
302
        float fip0            = -xr*(prr3*atomJ.q-prr5*dkr+prr7*qkr)
                             - prr3*atomJ.labFrameDipole_X + 2.0f*prr5*qkx;
303

Peter Eastman's avatar
Peter Eastman committed
304
305
        float fip1            = -yr*(prr3*atomJ.q-prr5*dkr+prr7*qkr)
                             - prr3*atomJ.labFrameDipole_Y + 2.0f*prr5*qky;
306

Peter Eastman's avatar
Peter Eastman committed
307
308
        float fip2            = -zr*(prr3*atomJ.q-prr5*dkr+prr7*qkr)
                             - prr3*atomJ.labFrameDipole_Z + 2.0f*prr5*qkz;
309

Peter Eastman's avatar
Peter Eastman committed
310
311
        float fkp0            = xr*(prr3*atomI.q+prr5*dir+prr7*qir)
                             - prr3*atomI.labFrameDipole_X - 2.0f*prr5*qix;
312

Peter Eastman's avatar
Peter Eastman committed
313
314
        float fkp1            = yr*(prr3*atomI.q+prr5*dir+prr7*qir)
                             - prr3*atomI.labFrameDipole_Y - 2.0f*prr5*qiy;
315

Peter Eastman's avatar
Peter Eastman committed
316
317
        float fkp2            = zr*(prr3*atomI.q+prr5*dir+prr7*qir)
                             - prr3*atomI.labFrameDipole_Z - 2.0f*prr5*qiz;
318

Peter Eastman's avatar
Peter Eastman committed
319
        // increment the field at each site due to this interaction
320

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

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

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

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

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

355
356
357
358
359
360
361
362
363
364
#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
365
366
367
368
369
370

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

371
372
373
374
375
376
377
378
379
380
381
382
#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
383
384
385
386
387
388
389
390
391
/**---------------------------------------------------------------------------------------

   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
392
#ifdef AMOEBA_DEBUG
Mark Friedrichs's avatar
Mark Friedrichs committed
393
394
395
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
396
#endif
Mark Friedrichs's avatar
Mark Friedrichs committed
397

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

   Compute fixed electric field using PME

   @param amoebaGpu        amoebaGpu context

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

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

#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
422
    memset( debugArray->_pSysStream[0], 0, sizeof( float )*4*paddedNumberOfAtoms*paddedNumberOfAtoms);
423
424
425
426
    debugArray->Upload();

    // print intermediate results for the targetAtom 

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

430
431
    kClearFields_3( amoebaGpu, 2 );

Mark Friedrichs's avatar
Mark Friedrichs committed
432
433
434
435
436
437
438
    // 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)
439
            maxThreads = 192;
Mark Friedrichs's avatar
Mark Friedrichs committed
440
441
442
443
444
        else
            maxThreads = 64;
        threadsPerBlock = std::min(getThreadsPerBlock(amoebaGpu, sizeof(FixedFieldParticle)), maxThreads);
    }    

445
    if (gpu->bOutputBufferPerWarp){
Mark Friedrichs's avatar
Mark Friedrichs committed
446
447
        kCalculateAmoebaPmeDirectFixedE_FieldN2ByWarp_kernel<<<amoebaGpu->nonbondBlocks, threadsPerBlock, sizeof(FixedFieldParticle)*threadsPerBlock>>>(
                                                                           gpu->sim.pInteractingWorkUnit,
448
449
450
451
452
453
454
455
                                                                           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
456
457
        kCalculateAmoebaPmeDirectFixedE_FieldN2_kernel<<<amoebaGpu->nonbondBlocks, threadsPerBlock, sizeof(FixedFieldParticle)*threadsPerBlock>>>(
                                                                           gpu->sim.pInteractingWorkUnit,
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
                                                                           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
473
474
475
        (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
476
        (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
477
                        amoebaGpu->nonbondBlocks, threadsPerBlock, amoebaGpu->bOutputBufferPerWarp,
Mark Friedrichs's avatar
Mark Friedrichs committed
478
479
                        sizeof(FixedFieldParticle), sizeof(FixedFieldParticle)*amoebaGpu->nonbondThreadsPerBlock, amoebaGpu->energyOutputBuffers, 
                        (*gpu->psInteractionCount)[0], gpu->sim.workUnits, gpu->bOutputBufferPerWarp );
480
        (void) fflush( amoebaGpu->log );
Mark Friedrichs's avatar
Mark Friedrichs committed
481
482
/*
        (void) fprintf( amoebaGpu->log, "Out WorkArray_3_[1,2]  paddedNumberOfAtoms=%d\n",  amoebaGpu->paddedNumberOfAtoms, amoebaGpu->outputBuffers );
483
484
        amoebaGpu->psWorkArray_3_1->Download();
        amoebaGpu->psWorkArray_3_2->Download();
Mark Friedrichs's avatar
Mark Friedrichs committed
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
        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
510
511
512
*/
        amoebaGpu->psE_Field->Download();
        amoebaGpu->psE_FieldPolar->Download();
Mark Friedrichs's avatar
Mark Friedrichs committed
513
514
        (void) fprintf( amoebaGpu->log,"E-field (includes self term)" );
        int maxPrint             = 3002;
515
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
        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
546
        amoebaGpu->gpuContext->psPosq4->Download();
547
548
        for( int jj = 0; jj < gpu->natoms; jj++ ){
            int debugIndex = jj;
Mark Friedrichs's avatar
Mark Friedrichs committed
549
if( fabs(debugArray->_pSysStream[0][jj+3*paddedNumberOfAtoms].x) > 0.0 ){
550
            (void) fprintf( amoebaGpu->log,"%5d PmeFixedEField\n", jj );
Mark Friedrichs's avatar
Mark Friedrichs committed
551
            for( int kk = 0; kk < 7; kk++ ){
552
553
554
555
556
557
                (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
558
}
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576

        }

        // 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
577
        if( 0 ){
578
579
580
581
582
583
584
585
586
587
            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
588
589
590

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