kCalculateAmoebaCudaRotateFrame.cu 22.7 KB
Newer Older
Mark Friedrichs's avatar
Mark Friedrichs committed
1
2
3
4
//-----------------------------------------------------------------------------------------

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

Mark Friedrichs's avatar
Mark Friedrichs committed
5
#include "cudaKernels.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
6
#include "amoebaCudaKernels.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
7
#include "kCalculateAmoebaCudaUtilities.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
8
9
10
11
12
13
14
15
16
17

#include <stdio.h>
#include <cuda.h>
#include <cstdlib>
using namespace std; 

#define SQRT sqrtf

static __constant__ cudaGmxSimulation cSim;
static __constant__ cudaAmoebaGmxSimulation cAmoebaSim;
18
extern __global__ void kFindInteractionsWithinBlocksPeriodic_kernel(unsigned int*);
Mark Friedrichs's avatar
Mark Friedrichs committed
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

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

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

40
__device__ static float normVector3( float* vector )
Mark Friedrichs's avatar
Mark Friedrichs committed
41
42
43
44
45
46
47
48
49
50
51
52
53
{

    float norm                    = DOT3( vector, vector );
    float returnNorm              = SQRT( norm );
    norm                          = returnNorm > 0.0f ? 1.0f/returnNorm : 0.0f;

    vector[0]                    *= norm;
    vector[1]                    *= norm;
    vector[2]                    *= norm;

    return returnNorm;
}

Mark Friedrichs's avatar
Mark Friedrichs committed
54
55
#undef AMOEBA_DEBUG

56
57
58
59
60
61
// ZThenX     == 0
// Bisector   == 1
// ZBisect    == 2
// ThreeFold  == 3
// ZOnly      == 4
// NoAxisType == 5
62

Mark Friedrichs's avatar
Mark Friedrichs committed
63
64
65
__global__
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
66
#elif (__CUDA_ARCH__ >= 120)
Mark Friedrichs's avatar
Mark Friedrichs committed
67
68
69
70
__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
#endif
71
void kCudaComputeCheckChiral_kernel( void )
Mark Friedrichs's avatar
Mark Friedrichs committed
72
73
{

74
75
76
77
78
79
    const int AD          = 0;
    const int BD          = 1;
    const int CD          = 2;
    const int C           = 3;
    float delta[4][3];
 
80
81
    float4* particleCoord        = cSim.pPosq;
    int4* multiPoleParticles     = cAmoebaSim.pMultipoleParticlesIdsAndAxisType;
82
83
84
85
86
    float* labFrameDipole        = cAmoebaSim.pLabFrameDipole;
    float* labFrameQuadrupole    = cAmoebaSim.pLabFrameQuadrupole;
 
    // ---------------------------------------------------------------------------------------
 
Mark Friedrichs's avatar
Mark Friedrichs committed
87
    int particleIndex            = blockIdx.x*blockDim.x + threadIdx.x;
88
89
90
91
92
    int numberOfParticles        = cSim.atoms;
    while( particleIndex < numberOfParticles )
    { 
        // skip z-then-x
    
Mark Friedrichs's avatar
Mark Friedrichs committed
93
94
        int axisType             = multiPoleParticles[particleIndex].w; 
        if( axisType != 0 && multiPoleParticles[particleIndex].x >= 0 && multiPoleParticles[particleIndex].y >=0 && multiPoleParticles[particleIndex].z >= 0 )
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
131
        {
     
            // ---------------------------------------------------------------------------------------
         
            int particleA                = particleIndex;
            int particleB                = multiPoleParticles[particleIndex].z;
            int particleC                = multiPoleParticles[particleIndex].x;
            int particleD                = multiPoleParticles[particleIndex].y;
        
            delta[AD][0]                 = particleCoord[particleA].x - particleCoord[particleD].x;
            delta[AD][1]                 = particleCoord[particleA].y - particleCoord[particleD].y;
            delta[AD][2]                 = particleCoord[particleA].z - particleCoord[particleD].z;
        
            delta[BD][0]                 = particleCoord[particleB].x - particleCoord[particleD].x;
            delta[BD][1]                 = particleCoord[particleB].y - particleCoord[particleD].y;
            delta[BD][2]                 = particleCoord[particleB].z - particleCoord[particleD].z;
        
            delta[CD][0]                 = particleCoord[particleC].x - particleCoord[particleD].x;
            delta[CD][1]                 = particleCoord[particleC].y - particleCoord[particleD].y;
            delta[CD][2]                 = particleCoord[particleC].z - particleCoord[particleD].z;
        
            delta[C][0]                  = delta[BD][1]*delta[CD][2] - delta[BD][2]*delta[CD][1];
            delta[C][1]                  = delta[CD][1]*delta[AD][2] - delta[CD][2]*delta[AD][1];
            delta[C][2]                  = delta[AD][1]*delta[BD][2] - delta[AD][2]*delta[BD][1];
         
            float volume                 = delta[C][0]*delta[AD][0] + delta[C][1]*delta[BD][0] + delta[C][2]*delta[CD][0];
            if( volume < 0.0 ){
                labFrameDipole[particleIndex*3+1]            *= -1.0f; // pole(3,i)
                labFrameQuadrupole[particleIndex*9+1]        *= -1.0f; // pole(6,i)  && pole(8,i)
                labFrameQuadrupole[particleIndex*9+3]        *= -1.0f; // pole(10,i) && pole(12,i)
                labFrameQuadrupole[particleIndex*9+5]        *= -1.0f; // pole(6,i)  && pole(8,i)
                labFrameQuadrupole[particleIndex*9+7]        *= -1.0f; // pole(10,i) && pole(12,i)
            }
        }
    
        particleIndex                += gridDim.x*blockDim.x;
    }    
132
}
Mark Friedrichs's avatar
Mark Friedrichs committed
133

134
135
136
__global__
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
137
#elif (__CUDA_ARCH__ >= 120)
138
139
140
141
142
143
__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
#endif
void kCudaComputeLabFrameMoments_kernel( void )
{
Mark Friedrichs's avatar
Mark Friedrichs committed
144

145
146
147
148
    float vectorX[3];
    float vectorY[3];
    float vectorZ[3];
 
Mark Friedrichs's avatar
Mark Friedrichs committed
149
    int particleIndex            = blockIdx.x*blockDim.x + threadIdx.x;
150
151
152

    float4* particleCoord        = cSim.pPosq;
    int4* multiPoleParticles     = cAmoebaSim.pMultipoleParticlesIdsAndAxisType;
153
154
155
156
157
158
159
160
161
162
163
    float* labFrameDipole        = cAmoebaSim.pLabFrameDipole;
    float* labFrameQuadrupole    = cAmoebaSim.pLabFrameQuadrupole;
 
    // get coordinates of this atom and the z & x axis atoms
    // compute the vector between the atoms and 1/sqrt(d2), d2 is distance between
    // this atom and the axis atom
 
    // this atom is referred to as the k-atom in notes below
 
    // code common to ZThenX and Bisector
    
164
    while( particleIndex < cSim.atoms )
165
    {
Mark Friedrichs's avatar
Mark Friedrichs committed
166
167
168
169

        if( multiPoleParticles[particleIndex].x >= 0 && multiPoleParticles[particleIndex].z >= 0 )
        {
            float4 coordinatesThisParticle   = particleCoord[particleIndex];
170
         
Mark Friedrichs's avatar
Mark Friedrichs committed
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
            int multipoleParticleIndex       = multiPoleParticles[particleIndex].z;
            float4 coordinatesAxisParticle   = particleCoord[multipoleParticleIndex];
         
            vectorZ[0]                       = coordinatesAxisParticle.x - coordinatesThisParticle.x;
            vectorZ[1]                       = coordinatesAxisParticle.y - coordinatesThisParticle.y;
            vectorZ[2]                       = coordinatesAxisParticle.z - coordinatesThisParticle.z;
              
            multipoleParticleIndex           = multiPoleParticles[particleIndex].x; 
            coordinatesAxisParticle          = particleCoord[multipoleParticleIndex];
         
            vectorX[0]                       = coordinatesAxisParticle.x - coordinatesThisParticle.x;
            vectorX[1]                       = coordinatesAxisParticle.y - coordinatesThisParticle.y;
            vectorX[2]                       = coordinatesAxisParticle.z - coordinatesThisParticle.z;
         
            int axisType                     = multiPoleParticles[particleIndex].w; 
186
    
Mark Friedrichs's avatar
Mark Friedrichs committed
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
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
            /*
                z-only
                   (1) norm z
                   (2) select random x
                   (3) x = x - (x.z)z
                   (4) norm x
        
                z-then-x
                   (1) norm z
                   (2) norm x (not needed)
                   (3) x = x - (x.z)z
                   (4) norm x
        
                bisector
                   (1) norm z
                   (2) norm x 
                   (3) z = x + z
                   (4) norm z
                   (5) x = x - (x.z)z 
                   (6) norm x 
        
                z-bisect
                   (1) norm z
                   (2) norm x 
                   (3) norm y 
                   (3) x = x + y
                   (4) norm x
                   (5) x = x - (x.z)z 
                   (6) norm x 
        
                3-fold
                   (1) norm z
                   (2) norm x 
                   (3) norm y 
                   (4) z = x + y + z
                   (5) norm z
                   (6) x = x - (x.z)z 
                   (7) norm x 
        
            */
        
            // branch based on axis type
             
            float sum                   = normVector3( vectorZ );
        
            if( axisType == 1 ){
        
                // bisector
                
                sum                     = normVector3( vectorX );
                
                vectorZ[0]             += vectorX[0];
                vectorZ[1]             += vectorX[1];
                vectorZ[2]             += vectorX[2];
           
                sum                     = normVector3( vectorZ );
        
            } else if( axisType == 2 || axisType == 3 ){ 
         
                // z-bisect
        
                multipoleParticleIndex  = multiPoleParticles[particleIndex].y; 
                if( multipoleParticleIndex >= 0 && multipoleParticleIndex < cSim.atoms ){
                    coordinatesAxisParticle = particleCoord[multipoleParticleIndex];
                    vectorY[0]              = coordinatesAxisParticle.x - coordinatesThisParticle.x;
                    vectorY[1]              = coordinatesAxisParticle.y - coordinatesThisParticle.y;
                    vectorY[2]              = coordinatesAxisParticle.z - coordinatesThisParticle.z;
254
            
Mark Friedrichs's avatar
Mark Friedrichs committed
255
256
                    sum                     = normVector3( vectorY );
                    sum                     = normVector3( vectorX );
257
            
Mark Friedrichs's avatar
Mark Friedrichs committed
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
                    if( axisType == 2 ){
            
                        vectorX[0]         += vectorY[0];
                        vectorX[1]         += vectorY[1];
                        vectorX[2]         += vectorY[2];
                        sum                 = normVector3( vectorX );
             
                    } else { 
             
                        // 3-fold
                
                        vectorZ[0]         += vectorX[0] + vectorY[0];
                        vectorZ[1]         += vectorX[1] + vectorY[1];
                        vectorZ[2]         += vectorX[2] + vectorY[2];
                        sum                 = normVector3( vectorZ );
                    }
                }
         
            } else if( axisType >= 4 ){ 
277
        
Mark Friedrichs's avatar
Mark Friedrichs committed
278
279
280
                vectorX[0]             = 0.1f;
                vectorX[1]             = 0.1f;
                vectorX[2]             = 0.1f;
281
            }
Mark Friedrichs's avatar
Mark Friedrichs committed
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
            
            // x = x - (x.z)z
        
            float dot         = vectorZ[0]*vectorX[0] + vectorZ[1]*vectorX[1] + vectorZ[2]*vectorX[2];
                
            vectorX[0]       -= dot*vectorZ[0];
            vectorX[1]       -= dot*vectorZ[1];
            vectorX[2]       -= dot*vectorZ[2];
             
            sum               = normVector3( vectorX );
        
            vectorY[0]        = (vectorZ[1]*vectorX[2]) - (vectorZ[2]*vectorX[1]);
            vectorY[1]        = (vectorZ[2]*vectorX[0]) - (vectorZ[0]*vectorX[2]);
            vectorY[2]        = (vectorZ[0]*vectorX[1]) - (vectorZ[1]*vectorX[0]);
         
            // use identity rotation matrix for unrecognized axis types
        
            if( axisType < 0 || axisType > 4 ){
300
        
Mark Friedrichs's avatar
Mark Friedrichs committed
301
302
303
304
305
306
307
308
309
310
311
312
313
314
                vectorX[0] = 1.0f;
                vectorX[1] = 0.0f;
                vectorX[2] = 0.0f;
        
                vectorY[0] = 0.0f;
                vectorY[1] = 1.0f;
                vectorY[2] = 0.0f;
        
                vectorZ[0] = 0.0f;
                vectorZ[1] = 0.0f;
                vectorZ[2] = 1.0f;
            }
        
            unsigned int offset           = 3*particleIndex;
315
    
Mark Friedrichs's avatar
Mark Friedrichs committed
316
317
318
319
            float molDipole[3];
            molDipole[0]                  = labFrameDipole[offset];
            molDipole[1]                  = labFrameDipole[offset+1];
            molDipole[2]                  = labFrameDipole[offset+2];
320
            
Mark Friedrichs's avatar
Mark Friedrichs committed
321
            // set out-of-range elements to 0.0f
322
         
Mark Friedrichs's avatar
Mark Friedrichs committed
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
            labFrameDipole[offset]        = molDipole[0]*vectorX[0] + molDipole[1]*vectorY[0] + molDipole[2]*vectorZ[0];
            labFrameDipole[offset+1]      = molDipole[0]*vectorX[1] + molDipole[1]*vectorY[1] + molDipole[2]*vectorZ[1];
            labFrameDipole[offset+2]      = molDipole[0]*vectorX[2] + molDipole[1]*vectorY[2] + molDipole[2]*vectorZ[2];
            
            // ---------------------------------------------------------------------------------------
            
            float mPole[3][3];
            offset                        = 9*particleIndex;
            
            mPole[0][0]                   = labFrameQuadrupole[offset];
            mPole[0][1]                   = labFrameQuadrupole[offset+1];
            mPole[0][2]                   = labFrameQuadrupole[offset+2];
        
            mPole[1][0]                   = labFrameQuadrupole[offset+3];
            mPole[1][1]                   = labFrameQuadrupole[offset+4];
            mPole[1][2]                   = labFrameQuadrupole[offset+5];
        
            mPole[2][0]                   = labFrameQuadrupole[offset+6];
            mPole[2][1]                   = labFrameQuadrupole[offset+7];
            mPole[2][2]                   = labFrameQuadrupole[offset+8];
        
            labFrameQuadrupole[offset+8]  = vectorX[2]*(vectorX[2]*mPole[0][0] + vectorY[2]*mPole[0][1] + vectorZ[2]*mPole[0][2]);
            labFrameQuadrupole[offset+8] += vectorY[2]*(vectorX[2]*mPole[1][0] + vectorY[2]*mPole[1][1] + vectorZ[2]*mPole[1][2]);
            labFrameQuadrupole[offset+8] += vectorZ[2]*(vectorX[2]*mPole[2][0] + vectorY[2]*mPole[2][1] + vectorZ[2]*mPole[2][2]);
347
    
Mark Friedrichs's avatar
Mark Friedrichs committed
348
349
350
            labFrameQuadrupole[offset+4]  = vectorX[1]*(vectorX[1]*mPole[0][0] + vectorY[1]*mPole[0][1] + vectorZ[1]*mPole[0][2]);
            labFrameQuadrupole[offset+4] += vectorY[1]*(vectorX[1]*mPole[1][0] + vectorY[1]*mPole[1][1] + vectorZ[1]*mPole[1][2]);
            labFrameQuadrupole[offset+4] += vectorZ[1]*(vectorX[1]*mPole[2][0] + vectorY[1]*mPole[2][1] + vectorZ[1]*mPole[2][2]);
351
    
Mark Friedrichs's avatar
Mark Friedrichs committed
352
353
354
            labFrameQuadrupole[offset+5]  = vectorX[1]*(vectorX[2]*mPole[0][0] + vectorY[2]*mPole[0][1] + vectorZ[2]*mPole[0][2]);
            labFrameQuadrupole[offset+5] += vectorY[1]*(vectorX[2]*mPole[1][0] + vectorY[2]*mPole[1][1] + vectorZ[2]*mPole[1][2]);
            labFrameQuadrupole[offset+5] += vectorZ[1]*(vectorX[2]*mPole[2][0] + vectorY[2]*mPole[2][1] + vectorZ[2]*mPole[2][2]);
355
    
Mark Friedrichs's avatar
Mark Friedrichs committed
356
357
358
            labFrameQuadrupole[offset]    = vectorX[0]*(vectorX[0]*mPole[0][0] + vectorY[0]*mPole[0][1] + vectorZ[0]*mPole[0][2]);
            labFrameQuadrupole[offset]   += vectorY[0]*(vectorX[0]*mPole[1][0] + vectorY[0]*mPole[1][1] + vectorZ[0]*mPole[1][2]);
            labFrameQuadrupole[offset]   += vectorZ[0]*(vectorX[0]*mPole[2][0] + vectorY[0]*mPole[2][1] + vectorZ[0]*mPole[2][2]);
359
    
Mark Friedrichs's avatar
Mark Friedrichs committed
360
361
362
            labFrameQuadrupole[offset+1]  = vectorX[0]*(vectorX[1]*mPole[0][0] + vectorY[1]*mPole[0][1] + vectorZ[1]*mPole[0][2]);
            labFrameQuadrupole[offset+1] += vectorY[0]*(vectorX[1]*mPole[1][0] + vectorY[1]*mPole[1][1] + vectorZ[1]*mPole[1][2]);
            labFrameQuadrupole[offset+1] += vectorZ[0]*(vectorX[1]*mPole[2][0] + vectorY[1]*mPole[2][1] + vectorZ[1]*mPole[2][2]);
363
    
Mark Friedrichs's avatar
Mark Friedrichs committed
364
365
366
            labFrameQuadrupole[offset+2]  = vectorX[0]*(vectorX[2]*mPole[0][0] + vectorY[2]*mPole[0][1] + vectorZ[2]*mPole[0][2]);
            labFrameQuadrupole[offset+2] += vectorY[0]*(vectorX[2]*mPole[1][0] + vectorY[2]*mPole[1][1] + vectorZ[2]*mPole[1][2]);
            labFrameQuadrupole[offset+2] += vectorZ[0]*(vectorX[2]*mPole[2][0] + vectorY[2]*mPole[2][1] + vectorZ[2]*mPole[2][2]);
Mark Friedrichs's avatar
Mark Friedrichs committed
367
     
Mark Friedrichs's avatar
Mark Friedrichs committed
368
369
370
            labFrameQuadrupole[offset+3]  = labFrameQuadrupole[offset+1];
            labFrameQuadrupole[offset+6]  = labFrameQuadrupole[offset+2];
            labFrameQuadrupole[offset+7]  = labFrameQuadrupole[offset+5];
371
    
Mark Friedrichs's avatar
Mark Friedrichs committed
372
        }
Mark Friedrichs's avatar
Mark Friedrichs committed
373

374
        particleIndex                += gridDim.x*blockDim.x;
375
    }
376
     
Mark Friedrichs's avatar
Mark Friedrichs committed
377
378
379
380
381
382
383
384
385
386
387
388
389
}

void cudaComputeAmoebaLabFrameMoments( amoebaGpuContext amoebaGpu )
{

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

   static const char* methodName = "computeCudaAmoebaLabFrameMoments";

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

    gpuContext gpu    = amoebaGpu->gpuContext;

390
    int numBlocks     = gpu->sim.blocks;
391
    int numThreads    = gpu->sim.threads_per_block;
Mark Friedrichs's avatar
Mark Friedrichs committed
392
393
394

//#define AMOEBA_DEBUG  
#ifdef AMOEBA_DEBUG
Mark Friedrichs's avatar
Mark Friedrichs committed
395
    if( amoebaGpu->log ){
Mark Friedrichs's avatar
Mark Friedrichs committed
396
397
398
        (void) fprintf( amoebaGpu->log, "%s: numBlocks/atoms=%d\n", methodName, numBlocks ); (void) fflush( amoebaGpu->log );
        amoebaGpu->psMultipoleParticlesIdsAndAxisType->Download();
        amoebaGpu->psMolecularDipole->Download();
399
        amoebaGpu->psMultipoleParticlesTorqueBufferIndices->Download();
Mark Friedrichs's avatar
Mark Friedrichs committed
400
401
402
        gpu->psPosq4->Download();
        for( int ii = 0; ii < gpu->natoms; ii++ ){
            int mIndex = 3*ii;
403
             (void) fprintf( amoebaGpu->log,"%6d [%6d %6d %6d %6d] x[%16.9e %16.9e %16.9e] %s [%6d %6d %6d %6d]\n", ii,
Mark Friedrichs's avatar
Mark Friedrichs committed
404
405
                             amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysData[ii].x,
                             amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysData[ii].y,
Mark Friedrichs's avatar
Mark Friedrichs committed
406
                             amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysData[ii].z,
Mark Friedrichs's avatar
Mark Friedrichs committed
407
408
409
                             amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysData[ii].w,
                             gpu->psPosq4->_pSysData[ii].x,
                             gpu->psPosq4->_pSysData[ii].y,
410
411
412
413
414
                             gpu->psPosq4->_pSysData[ii].z, (amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysData[ii].w > 1 ? " XXX" : ""),
                             amoebaGpu->psMultipoleParticlesTorqueBufferIndices->_pSysData[ii].x,
                             amoebaGpu->psMultipoleParticlesTorqueBufferIndices->_pSysData[ii].y,
                             amoebaGpu->psMultipoleParticlesTorqueBufferIndices->_pSysData[ii].z,
                             amoebaGpu->psMultipoleParticlesTorqueBufferIndices->_pSysData[ii].w );
Mark Friedrichs's avatar
Mark Friedrichs committed
415
            //if( ii == 30 )ii = gpu->natoms - 30;
Mark Friedrichs's avatar
Mark Friedrichs committed
416
417
418
419
        }
    }
#endif

420
421
422
423
424
425
426
    // copy molecular moments to lab frame moment arrays
    // check if chiral center requires moments to have sign flipped
    // compute lab frame moments

    cudaMemcpy( amoebaGpu->psLabFrameDipole->_pDevData,      amoebaGpu->psMolecularDipole->_pDevData,    3*gpu->sim.paddedNumberOfAtoms*sizeof( float ), cudaMemcpyDeviceToDevice ); 
    cudaMemcpy( amoebaGpu->psLabFrameQuadrupole->_pDevData, amoebaGpu->psMolecularQuadrupole->_pDevData, 9*gpu->sim.paddedNumberOfAtoms*sizeof( float ), cudaMemcpyDeviceToDevice ); 

427
428
429
430
    kCudaComputeCheckChiral_kernel<<< numBlocks, numThreads>>> ( );
    LAUNCHERROR("kCudaComputeCheckChiral");

    kCudaComputeLabFrameMoments_kernel<<< numBlocks, numThreads>>> ( );
Mark Friedrichs's avatar
Mark Friedrichs committed
431
432
433
434
435
436
437
438
439
440
441
442
    LAUNCHERROR(methodName);

}

void kCalculateAmoebaMultipoleForces(amoebaGpuContext amoebaGpu, bool hasAmoebaGeneralizedKirkwood ) 
{
    std::string methodName = "kCalculateAmoebaMultipoleForces";

    // compute lab frame moments

    cudaComputeAmoebaLabFrameMoments( amoebaGpu );

Mark Friedrichs's avatar
Mark Friedrichs committed
443
    if( 0 ){
Mark Friedrichs's avatar
Mark Friedrichs committed
444
445
446
447
        gpuContext gpu                       = amoebaGpu->gpuContext;
        std::vector<int> fileId;
        //fileId.push_back( 0 );
        VectorOfDoubleVectors outputVector;
448
449
450
        //cudaLoadCudaFloat4Array( gpu->natoms, 3, gpu->psPosq4,              outputVector, gpu->psAtomIndex->_pSysData, 1.0f );
        cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psLabFrameDipole,     outputVector, gpu->psAtomIndex->_pSysData, 1.0f );
        cudaLoadCudaFloatArray( gpu->natoms,  9, amoebaGpu->psLabFrameQuadrupole, outputVector, gpu->psAtomIndex->_pSysData, 1.0f );
Mark Friedrichs's avatar
Mark Friedrichs committed
451
452
453
        cudaWriteVectorOfDoubleVectorsToFile( "CudaLabMoments", fileId, outputVector );
    }   

Mark Friedrichs's avatar
Mark Friedrichs committed
454
455
456
457
458
459
    // compute fixed E-field and mutual induced field 

    if( hasAmoebaGeneralizedKirkwood ){
        cudaComputeAmoebaFixedEAndGkFields( amoebaGpu );
        cudaComputeAmoebaMutualInducedAndGkField( amoebaGpu );
    } else {
460

461
        if( amoebaGpu->multipoleNonbondedMethod == AMOEBA_NO_CUTOFF ){
462

463
464
            cudaComputeAmoebaFixedEField( amoebaGpu );
            cudaComputeAmoebaMutualInducedField( amoebaGpu );
465

466
        } else {
467

Mark Friedrichs's avatar
Mark Friedrichs committed
468
469
470
471
472
            gpuContext gpu = amoebaGpu->gpuContext;
            kFindBlockBoundsPeriodic_kernel<<<(gpu->psGridBoundingBox->_length+63)/64, 64>>>();
            LAUNCHERROR("kFindBlockBoundsPeriodic");
            kFindBlocksWithInteractionsPeriodic_kernel<<<gpu->sim.interaction_blocks, gpu->sim.interaction_threads_per_block>>>();
            LAUNCHERROR("kFindBlocksWithInteractionsPeriodic");
Mark Friedrichs's avatar
Mark Friedrichs committed
473
            //compactStream(gpu->compactPlan, gpu->sim.pInteractingWorkUnit, gpu->sim.pWorkUnit, gpu->sim.pInteractionFlag, gpu->sim.workUnits, gpu->sim.pInteractionCount);
Mark Friedrichs's avatar
Mark Friedrichs committed
474
            compactStream(gpu->compactPlan, gpu->sim.pInteractingWorkUnit, amoebaGpu->psWorkUnit->_pDevData, gpu->sim.pInteractionFlag, gpu->sim.workUnits, gpu->sim.pInteractionCount);
Mark Friedrichs's avatar
Mark Friedrichs committed
475
476
477
            kFindInteractionsWithinBlocksPeriodic_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
                    sizeof(unsigned int)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
            LAUNCHERROR("kFindInteractionsWithinBlocksPeriodic");
478

479
            cudaComputeAmoebaPmeFixedEField( amoebaGpu );
480
            cudaComputeAmoebaPmeMutualInducedField( amoebaGpu );
481
        }
Mark Friedrichs's avatar
Mark Friedrichs committed
482
483
484
485
    }

    // check if induce dipole calculation converged -- abort if it did not

486
    if( amoebaGpu->mutualInducedDone == 0 ){
Mark Friedrichs's avatar
Mark Friedrichs committed
487
488
       (void) fprintf( stderr, "%s induced dipole calculation did not converge -- aborting!\n", methodName.c_str() );
       (void) fflush( stderr );
Mark Friedrichs's avatar
Mark Friedrichs committed
489
490
491
492
493
       exit(-1);
    }

    // calculate electrostatic forces

494
    if( amoebaGpu->multipoleNonbondedMethod == AMOEBA_NO_CUTOFF ){
495
        cudaComputeAmoebaElectrostatic( amoebaGpu, (hasAmoebaGeneralizedKirkwood ? 0 : 1) );
496
    } else {
Mark Friedrichs's avatar
Mark Friedrichs committed
497
        cudaComputeAmoebaPmeElectrostatic( amoebaGpu );
498
    }
499

Mark Friedrichs's avatar
Mark Friedrichs committed
500
501
502
}

#undef AMOEBA_DEBUG