kCalculateAmoebaCudaPmeDirectElectrostatic.cu 50.7 KB
Newer Older
Mark Friedrichs's avatar
Mark Friedrichs committed
1
///-----------------------------------------------------------------------------------------
2
3
4
5

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

#include "amoebaGpuTypes.h"
6
#include "cudaKernels.h"
7
8
9
#include "amoebaCudaKernels.h"
#include "kCalculateAmoebaCudaUtilities.h"

10
//#define AMOEBA_DEBUG
11
12
13
14

static __constant__ cudaGmxSimulation cSim;
static __constant__ cudaAmoebaGmxSimulation cAmoebaSim;

Mark Friedrichs's avatar
Mark Friedrichs committed
15
void SetCalculateAmoebaPmeDirectElectrostaticSim(amoebaGpuContext amoebaGpu)
16
17
18
{
    cudaError_t status;
    gpuContext gpu = amoebaGpu->gpuContext;
Peter Eastman's avatar
Peter Eastman committed
19
    status         = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));
Mark Friedrichs's avatar
Mark Friedrichs committed
20
    RTERROR(status, "SetCalculateAmoebaPmeDirectElectrostaticSim: cudaMemcpyToSymbol: SetSim copy to cSim failed");
Peter Eastman's avatar
Peter Eastman committed
21
    status         = cudaMemcpyToSymbol(cAmoebaSim, &amoebaGpu->amoebaSim, sizeof(cudaAmoebaGmxSimulation));
Mark Friedrichs's avatar
Mark Friedrichs committed
22
    RTERROR(status, "SetCalculateAmoebaPmeDirectElectrostaticSim: cudaMemcpyToSymbol: SetSim copy to cAmoebaSim failed");
23
24
}

Mark Friedrichs's avatar
Mark Friedrichs committed
25
void GetCalculateAmoebaPmeDirectElectrostaticSim(amoebaGpuContext amoebaGpu)
26
27
28
{
    cudaError_t status;
    gpuContext gpu = amoebaGpu->gpuContext;
Peter Eastman's avatar
Peter Eastman committed
29
    status = cudaMemcpyFromSymbol(&gpu->sim, cSim, sizeof(cudaGmxSimulation));
Mark Friedrichs's avatar
Mark Friedrichs committed
30
    RTERROR(status, "GetCalculateAmoebaPmeDirectElectrostaticSim: cudaMemcpyFromSymbol: SetSim copy from cSim failed");
Peter Eastman's avatar
Peter Eastman committed
31
    status = cudaMemcpyFromSymbol(&amoebaGpu->amoebaSim, cAmoebaSim, sizeof(cudaAmoebaGmxSimulation));
Mark Friedrichs's avatar
Mark Friedrichs committed
32
    RTERROR(status, "GetCalculateAmoebaPmeDirectElectrostaticSim: cudaMemcpyFromSymbol: SetSim copy from cAmoebaSim failed");
33
34
}

Peter Eastman's avatar
Peter Eastman committed
35
36
37
static int const PScaleIndex            =  0;
static int const DScaleIndex            =  1;
static int const UScaleIndex            =  2;
38
static int const MScaleIndex            =  3;
Mark Friedrichs's avatar
Mark Friedrichs committed
39
static int const LastScalingIndex       =  4;
40

Mark Friedrichs's avatar
Mark Friedrichs committed
41
struct PmeDirectElectrostaticParticle {
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

    // coordinates charge

    float x;
    float y;
    float z;
    float q;

    // lab frame dipole

    float labFrameDipole[3];

    // lab frame quadrupole

    float labFrameQuadrupole[9];

    // induced dipole

    float inducedDipole[3];

    // polar induced dipole

    float inducedDipoleP[3];

    // scaling factors

    float thole;
    float damp;

    float force[3];

    float torque[3];
    float padding;

Mark Friedrichs's avatar
Mark Friedrichs committed
76
77
    float tempForce[3];
    float tempTorque[3];
78
79
};

Mark Friedrichs's avatar
Mark Friedrichs committed
80
81
82
83
84
85
86
87
88
89
__device__ void sumTempBuffer( PmeDirectElectrostaticParticle& atomI, PmeDirectElectrostaticParticle& atomJ ){

    atomI.tempForce[0]  += atomJ.tempForce[0];
    atomI.tempForce[1]  += atomJ.tempForce[1];
    atomI.tempForce[2]  += atomJ.tempForce[2];

    atomI.tempTorque[0] += atomJ.tempTorque[0];
    atomI.tempTorque[1] += atomJ.tempTorque[1];
    atomI.tempTorque[2] += atomJ.tempTorque[2];
}
Mark Friedrichs's avatar
Mark Friedrichs committed
90

91
92
93
94
/*
__device__ static void debugSetup( unsigned int atomI, unsigned int atomJ,
                                   float4* debugArray, float4* pullBack )
{
95
    unsigned int index                 = atomI + atomJ*cSim.paddedNumberOfAtoms;
96
97
98
99
100
101
102
103
    float blockId                      = 111.0f;

    debugArray[index].x                = (float) atomI;
    debugArray[index].y                = (float) atomJ;
    debugArray[index].z                = 0.0f;
    debugArray[index].w                = blockId;

    for( int pullIndex = 0; pullIndex < 1; pullIndex++ ){
104
        index                             += cSim.paddedNumberOfAtoms;
105
106
107
108
        debugArray[index].x                = pullBack[pullIndex].x;
        debugArray[index].y                = pullBack[pullIndex].y;
        debugArray[index].z                = pullBack[pullIndex].z;
        debugArray[index].w                = pullBack[pullIndex].w;
Peter Eastman's avatar
Peter Eastman committed
109
    }
110
111
112
113
}
*/

// self-energy for PME
Mark Friedrichs's avatar
Mark Friedrichs committed
114

Mark Friedrichs's avatar
Mark Friedrichs committed
115
116
117
__device__ static void calculatePmeSelfEnergyElectrostaticPairIxn_kernel( PmeDirectElectrostaticParticle& atomI, float* energy)
{
    float term   = 2.0f*cSim.alphaEwald*cSim.alphaEwald;
Mark Friedrichs's avatar
Mark Friedrichs committed
118
    float fterm  = -cSim.alphaEwald/cAmoebaSim.sqrtPi;
Mark Friedrichs's avatar
Mark Friedrichs committed
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

    float cii    = atomI.q*atomI.q;

    float dii    = atomI.labFrameDipole[0]*atomI.labFrameDipole[0] +
                   atomI.labFrameDipole[1]*atomI.labFrameDipole[1] +
                   atomI.labFrameDipole[2]*atomI.labFrameDipole[2];

    float qii    = atomI.labFrameQuadrupole[0]*atomI.labFrameQuadrupole[0] +
                   atomI.labFrameQuadrupole[4]*atomI.labFrameQuadrupole[4] +
                   atomI.labFrameQuadrupole[8]*atomI.labFrameQuadrupole[8] + 2.0f*(
                   atomI.labFrameQuadrupole[1]*atomI.labFrameQuadrupole[1] +
                   atomI.labFrameQuadrupole[2]*atomI.labFrameQuadrupole[2] +
                   atomI.labFrameQuadrupole[5]*atomI.labFrameQuadrupole[5]);

    float uii    = atomI.labFrameDipole[0]*atomI.inducedDipole[0] + atomI.labFrameDipole[1]*atomI.inducedDipole[1] + atomI.labFrameDipole[2]*atomI.inducedDipole[2];

    *energy      = (cii + term*(dii/3.0f + 2.0f*term*qii/5.0f));
    *energy     += term*uii/3.0f;
    *energy     *= fterm;
Mark Friedrichs's avatar
Mark Friedrichs committed
138
}
Mark Friedrichs's avatar
Mark Friedrichs committed
139
140
141
142
143

// self-torque for PME

__device__ static void calculatePmeSelfTorqueElectrostaticPairIxn_kernel( PmeDirectElectrostaticParticle& atomI)
{
Mark Friedrichs's avatar
Mark Friedrichs committed
144
    float term       = (4.0f/3.0f)*(cSim.alphaEwald*cSim.alphaEwald*cSim.alphaEwald)/cAmoebaSim.sqrtPi;
Mark Friedrichs's avatar
Mark Friedrichs committed
145
146
147
148
149

    float uix        = 0.5f*(atomI.inducedDipole[0] + atomI.inducedDipoleP[0]);
    float uiy        = 0.5f*(atomI.inducedDipole[1] + atomI.inducedDipoleP[1]);
    float uiz        = 0.5f*(atomI.inducedDipole[2] + atomI.inducedDipoleP[2]);

Mark Friedrichs's avatar
Mark Friedrichs committed
150
151
152
    atomI.torque[0] += term*(atomI.labFrameDipole[1]*uiz - atomI.labFrameDipole[2]*uiy);
    atomI.torque[1] += term*(atomI.labFrameDipole[2]*uix - atomI.labFrameDipole[0]*uiz);
    atomI.torque[2] += term*(atomI.labFrameDipole[0]*uiy - atomI.labFrameDipole[1]*uix);
Mark Friedrichs's avatar
Mark Friedrichs committed
153
154
155
}

__device__ void calculatePmeDirectElectrostaticPairIxn_kernel( PmeDirectElectrostaticParticle& atomI,   PmeDirectElectrostaticParticle& atomJ,
Mark Friedrichs's avatar
Mark Friedrichs committed
156
                                                               float* scalingFactors, float4 forceTorqueEnergy[3]
157
#ifdef AMOEBA_DEBUG
Peter Eastman's avatar
Peter Eastman committed
158
                                                               ,float4* debugArray
159
160
#endif
 ){
Mark Friedrichs's avatar
Mark Friedrichs committed
161

162
163
164
165
    float xr    = atomJ.x - atomI.x;
    float yr    = atomJ.y - atomI.y;
    float zr    = atomJ.z - atomI.z;

Peter Eastman's avatar
Peter Eastman committed
166
    // periodic box
167
168
169
170
171
172

    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;

    float r2    = xr*xr + yr*yr + zr*zr;
Mark Friedrichs's avatar
Mark Friedrichs committed
173
    if( r2 <= cSim.nonbondedCutoffSqr ){
Mark Friedrichs's avatar
Mark Friedrichs committed
174

175
176
        float r      = sqrt(r2);
        float ck     = atomJ.q;
Peter Eastman's avatar
Peter Eastman committed
177
178
179
180
181

        // set the permanent multipole and induced dipole values;

        float ci    = atomI.q;

182
183
184
        float di1   = atomI.labFrameDipole[0];
        float di2   = atomI.labFrameDipole[1];
        float di3   = atomI.labFrameDipole[2];
Peter Eastman's avatar
Peter Eastman committed
185

186
187
188
189
190
191
192
193
194
        float qi1   = atomI.labFrameQuadrupole[0];
        float qi2   = atomI.labFrameQuadrupole[1];
        float qi3   = atomI.labFrameQuadrupole[2];
        float qi4   = atomI.labFrameQuadrupole[3];
        float qi5   = atomI.labFrameQuadrupole[4];
        float qi6   = atomI.labFrameQuadrupole[5];
        float qi7   = atomI.labFrameQuadrupole[6];
        float qi8   = atomI.labFrameQuadrupole[7];
        float qi9   = atomI.labFrameQuadrupole[8];
Peter Eastman's avatar
Peter Eastman committed
195

196
197
198
        float dk1  = atomJ.labFrameDipole[0];
        float dk2  = atomJ.labFrameDipole[1];
        float dk3  = atomJ.labFrameDipole[2];
Peter Eastman's avatar
Peter Eastman committed
199

200
201
202
203
204
205
206
207
208
        float qk1  = atomJ.labFrameQuadrupole[0];
        float qk2  = atomJ.labFrameQuadrupole[1];
        float qk3  = atomJ.labFrameQuadrupole[2];
        float qk4  = atomJ.labFrameQuadrupole[3];
        float qk5  = atomJ.labFrameQuadrupole[4];
        float qk6  = atomJ.labFrameQuadrupole[5];
        float qk7  = atomJ.labFrameQuadrupole[6];
        float qk8  = atomJ.labFrameQuadrupole[7];
        float qk9  = atomJ.labFrameQuadrupole[8];
Peter Eastman's avatar
Peter Eastman committed
209

Mark Friedrichs's avatar
Mark Friedrichs committed
210
        // calculate the real space error function terms
211

212
        float ralpha = cSim.alphaEwald*r;
Mark Friedrichs's avatar
Mark Friedrichs committed
213
        float bn0    = erfc(ralpha)/r;
214

215
        float alsq2  = 2.0f*cSim.alphaEwald*cSim.alphaEwald;
216
        float alsq2n = 0.0f;
Mark Friedrichs's avatar
Mark Friedrichs committed
217
        if( cSim.alphaEwald > 0.0f ){
218
            alsq2n = 1.0f/(cAmoebaSim.sqrtPi*cSim.alphaEwald);
219
220
221
222
        }
        float exp2a  = exp(-(ralpha*ralpha));

        alsq2n      *= alsq2;
223
        float bn1    = (bn0+alsq2n*exp2a)/r2;
224
225

        alsq2n      *= alsq2;
226
        float bn2    = (3.0f*bn1+alsq2n*exp2a)/r2;
227
228

        alsq2n      *= alsq2;
229
        float bn3    = (5.0f*bn2+alsq2n*exp2a)/r2;
230
231

        alsq2n      *= alsq2;
232
        float bn4    = (7.0f*bn3+alsq2n*exp2a)/r2;
233
234

        alsq2n      *= alsq2;
235
        float bn5    = (9.0f*bn4+alsq2n*exp2a)/r2;
236

Mark Friedrichs's avatar
Mark Friedrichs committed
237
        // apply Thole polarization damping to scale factors
238
239
240
241
242
243
244
245
246
247
248

        float rr1    = 1.0f/r;
        float rr3    = rr1 / r2;
        float rr5    = 3.0f * rr3 / r2;
        float rr7    = 5.0f * rr5 / r2;
        float rr9    = 7.0f * rr7 / r2;
        float rr11   = 9.0f * rr9 / r2;
        float scale3 = 1.0f;
        float scale5 = 1.0f;
        float scale7 = 1.0f;

249
250
251
        float ddsc31 = 0.0f;
        float ddsc32 = 0.0f;
        float ddsc33 = 0.0f;
252

253
254
255
        float ddsc51 = 0.0f;
        float ddsc52 = 0.0f;
        float ddsc53 = 0.0f;
256

257
258
259
        float ddsc71 = 0.0f;
        float ddsc72 = 0.0f;
        float ddsc73 = 0.0f;
260

261
        float damp   = atomI.damp*atomJ.damp;
262
        if( damp != 0.0f ){
263
            float pgamma = atomI.thole < atomJ.thole ? atomI.thole : atomJ.thole;
264
            float ratio  = r/damp;
Mark Friedrichs's avatar
Mark Friedrichs committed
265
                damp     = -pgamma*ratio*ratio*ratio;
266
267
268
269
270
271
272
273
            if( damp > -50.0f ){
                float expdamp  = exp(damp);
                scale3         = 1.0f - expdamp;
                   scale5      = 1.0f - (1.0f-damp)*expdamp;
                   scale7      = 1.0f - (1.0f-damp+0.6f*damp*damp)*expdamp;
                float temp3    = -3.0f * damp * expdamp / r2;
                float temp5    = -damp;
                float temp7    = -0.2f - 0.6f*damp;
Mark Friedrichs's avatar
Mark Friedrichs committed
274

Mark Friedrichs's avatar
Mark Friedrichs committed
275
276
277
                ddsc31         = temp3 * xr;
                ddsc32         = temp3 * yr;
                ddsc33         = temp3 * zr;
Mark Friedrichs's avatar
Mark Friedrichs committed
278

Mark Friedrichs's avatar
Mark Friedrichs committed
279
280
281
                ddsc51         = temp5 * ddsc31;
                ddsc52         = temp5 * ddsc32;
                ddsc53         = temp5 * ddsc33;
Mark Friedrichs's avatar
Mark Friedrichs committed
282

Mark Friedrichs's avatar
Mark Friedrichs committed
283
284
285
                ddsc71         = temp7 * ddsc51;
                ddsc72         = temp7 * ddsc52;
                ddsc73         = temp7 * ddsc53;
Mark Friedrichs's avatar
Mark Friedrichs committed
286
            }
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
        }

        float dsc3 = 1.0f - scale3*scalingFactors[DScaleIndex];
        float dsc5 = 1.0f - scale5*scalingFactors[DScaleIndex];
        float dsc7 = 1.0f - scale7*scalingFactors[DScaleIndex];

        float psc3 = 1.0f - scale3*scalingFactors[PScaleIndex];
        float psc5 = 1.0f - scale5*scalingFactors[PScaleIndex];
        float psc7 = 1.0f - scale7*scalingFactors[PScaleIndex];

        float usc3 = 1.0f - scale3*scalingFactors[UScaleIndex];
        float usc5 = 1.0f - scale5*scalingFactors[UScaleIndex];

        // construct necessary auxiliary vectors

302
303
304
305
        float dixdk1       = di2*dk3 - di3*dk2;
        float dixdk2       = di3*dk1 - di1*dk3;
        float dixdk3       = di1*dk2 - di2*dk1;

306
307
308
309
310
311
312
        float dixuk1       = di2*atomJ.inducedDipole[2]  - di3*atomJ.inducedDipole[1];
        float dixuk2       = di3*atomJ.inducedDipole[0]  - di1*atomJ.inducedDipole[2];
        float dixuk3       = di1*atomJ.inducedDipole[1]  - di2*atomJ.inducedDipole[0];
        float dkxui1       = dk2*atomI.inducedDipole[2]  - dk3*atomI.inducedDipole[1];
        float dkxui2       = dk3*atomI.inducedDipole[0]  - dk1*atomI.inducedDipole[2];
        float dkxui3       = dk1*atomI.inducedDipole[1]  - dk2*atomI.inducedDipole[0];

313
314
315
316
317
318
        float dixukp1      = di2*atomJ.inducedDipoleP[2] - di3*atomJ.inducedDipoleP[1];
        float dixukp2      = di3*atomJ.inducedDipoleP[0] - di1*atomJ.inducedDipoleP[2];
        float dixukp3      = di1*atomJ.inducedDipoleP[1] - di2*atomJ.inducedDipoleP[0];
        float dkxuip1      = dk2*atomI.inducedDipoleP[2] - dk3*atomI.inducedDipoleP[1];
        float dkxuip2      = dk3*atomI.inducedDipoleP[0] - dk1*atomI.inducedDipoleP[2];
        float dkxuip3      = dk1*atomI.inducedDipoleP[1] - dk2*atomI.inducedDipoleP[0];
319

320
321
322
        float dixr1        = di2*zr - di3*yr;
        float dixr2        = di3*xr - di1*zr;
        float dixr3        = di1*yr - di2*xr;
323

324
325
326
        float dkxr1        = dk2*zr - dk3*yr;
        float dkxr2        = dk3*xr - dk1*zr;
        float dkxr3        = dk1*yr - dk2*xr;
327

328
329
330
        float qir1         = qi1*xr + qi4*yr + qi7*zr;
        float qir2         = qi2*xr + qi5*yr + qi8*zr;
        float qir3         = qi3*xr + qi6*yr + qi9*zr;
331

332
333
334
        float qkr1         = qk1*xr + qk4*yr + qk7*zr;
        float qkr2         = qk2*xr + qk5*yr + qk8*zr;
        float qkr3         = qk3*xr + qk6*yr + qk9*zr;
335

336
337
338
        float qiqkr1       = qi1*qkr1 + qi4*qkr2 + qi7*qkr3;
        float qiqkr2       = qi2*qkr1 + qi5*qkr2 + qi8*qkr3;
        float qiqkr3       = qi3*qkr1 + qi6*qkr2 + qi9*qkr3;
339

340
341
342
        float qkqir1       = qk1*qir1 + qk4*qir2 + qk7*qir3;
        float qkqir2       = qk2*qir1 + qk5*qir2 + qk8*qir3;
        float qkqir3       = qk3*qir1 + qk6*qir2 + qk9*qir3;
343
344
345
346
347

        float qixqk1       = qi2*qk3 + qi5*qk6 + qi8*qk9 - qi3*qk2 - qi6*qk5 - qi9*qk8;
        float qixqk2       = qi3*qk1 + qi6*qk4 + qi9*qk7 - qi1*qk3 - qi4*qk6 - qi7*qk9;
        float qixqk3       = qi1*qk2 + qi4*qk5 + qi7*qk8 - qi2*qk1 - qi5*qk4 - qi8*qk7;

348
349
350
        float rxqir1       = yr*qir3 - zr*qir2;
        float rxqir2       = zr*qir1 - xr*qir3;
        float rxqir3       = xr*qir2 - yr*qir1;
351

352
353
354
        float rxqkr1       = yr*qkr3 - zr*qkr2;
        float rxqkr2       = zr*qkr1 - xr*qkr3;
        float rxqkr3       = xr*qkr2 - yr*qkr1;
355

356
357
358
        float rxqikr1      = yr*qiqkr3 - zr*qiqkr2;
        float rxqikr2      = zr*qiqkr1 - xr*qiqkr3;
        float rxqikr3      = xr*qiqkr2 - yr*qiqkr1;
359

360
361
362
        float rxqkir1      = yr*qkqir3 - zr*qkqir2;
        float rxqkir2      = zr*qkqir1 - xr*qkqir3;
        float rxqkir3      = xr*qkqir2 - yr*qkqir1;
363

364
365
366
        float qkrxqir1     = qkr2*qir3 - qkr3*qir2;
        float qkrxqir2     = qkr3*qir1 - qkr1*qir3;
        float qkrxqir3     = qkr1*qir2 - qkr2*qir1;
367

368
369
370
        float qidk1        = qi1*dk1 + qi4*dk2 + qi7*dk3;
        float qidk2        = qi2*dk1 + qi5*dk2 + qi8*dk3;
        float qidk3        = qi3*dk1 + qi6*dk2 + qi9*dk3;
371

372
373
374
        float qkdi1        = qk1*di1 + qk4*di2 + qk7*di3;
        float qkdi2        = qk2*di1 + qk5*di2 + qk8*di3;
        float qkdi3        = qk3*di1 + qk6*di2 + qk9*di3;
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391

        float qiuk1        = qi1*atomJ.inducedDipole[0]  + qi4*atomJ.inducedDipole[1]  + qi7*atomJ.inducedDipole[2];
        float qiuk2        = qi2*atomJ.inducedDipole[0]  + qi5*atomJ.inducedDipole[1]  + qi8*atomJ.inducedDipole[2];
        float qiuk3        = qi3*atomJ.inducedDipole[0]  + qi6*atomJ.inducedDipole[1]  + qi9*atomJ.inducedDipole[2];

        float qkui1        = qk1*atomI.inducedDipole[0]  + qk4*atomI.inducedDipole[1]  + qk7*atomI.inducedDipole[2];
        float qkui2        = qk2*atomI.inducedDipole[0]  + qk5*atomI.inducedDipole[1]  + qk8*atomI.inducedDipole[2];
        float qkui3        = qk3*atomI.inducedDipole[0]  + qk6*atomI.inducedDipole[1]  + qk9*atomI.inducedDipole[2];

        float qiukp1       = qi1*atomJ.inducedDipoleP[0] + qi4*atomJ.inducedDipoleP[1] + qi7*atomJ.inducedDipoleP[2];
        float qiukp2       = qi2*atomJ.inducedDipoleP[0] + qi5*atomJ.inducedDipoleP[1] + qi8*atomJ.inducedDipoleP[2];
        float qiukp3       = qi3*atomJ.inducedDipoleP[0] + qi6*atomJ.inducedDipoleP[1] + qi9*atomJ.inducedDipoleP[2];

        float qkuip1       = qk1*atomI.inducedDipoleP[0] + qk4*atomI.inducedDipoleP[1] + qk7*atomI.inducedDipoleP[2];
        float qkuip2       = qk2*atomI.inducedDipoleP[0] + qk5*atomI.inducedDipoleP[1] + qk8*atomI.inducedDipoleP[2];
        float qkuip3       = qk3*atomI.inducedDipoleP[0] + qk6*atomI.inducedDipoleP[1] + qk9*atomI.inducedDipoleP[2];

392
393
394
        float dixqkr1      = di2*qkr3 - di3*qkr2;
        float dixqkr2      = di3*qkr1 - di1*qkr3;
        float dixqkr3      = di1*qkr2 - di2*qkr1;
395

396
397
398
        float dkxqir1      = dk2*qir3 - dk3*qir2;
        float dkxqir2      = dk3*qir1 - dk1*qir3;
        float dkxqir3      = dk1*qir2 - dk2*qir1;
399

400
401
402
        float uixqkr1      = atomI.inducedDipole[1]*qkr3 - atomI.inducedDipole[2]*qkr2;
        float uixqkr2      = atomI.inducedDipole[2]*qkr1 - atomI.inducedDipole[0]*qkr3;
        float uixqkr3      = atomI.inducedDipole[0]*qkr2 - atomI.inducedDipole[1]*qkr1;
403

404
405
406
        float ukxqir1      = atomJ.inducedDipole[1]*qir3 - atomJ.inducedDipole[2]*qir2;
        float ukxqir2      = atomJ.inducedDipole[2]*qir1 - atomJ.inducedDipole[0]*qir3;
        float ukxqir3      = atomJ.inducedDipole[0]*qir2 - atomJ.inducedDipole[1]*qir1;
407

408
409
410
        float uixqkrp1     = atomI.inducedDipoleP[1]*qkr3 - atomI.inducedDipoleP[2]*qkr2;
        float uixqkrp2     = atomI.inducedDipoleP[2]*qkr1 - atomI.inducedDipoleP[0]*qkr3;
        float uixqkrp3     = atomI.inducedDipoleP[0]*qkr2 - atomI.inducedDipoleP[1]*qkr1;
411

412
413
414
        float ukxqirp1     = atomJ.inducedDipoleP[1]*qir3 - atomJ.inducedDipoleP[2]*qir2;
        float ukxqirp2     = atomJ.inducedDipoleP[2]*qir1 - atomJ.inducedDipoleP[0]*qir3;
        float ukxqirp3     = atomJ.inducedDipoleP[0]*qir2 - atomJ.inducedDipoleP[1]*qir1;
415

416
417
418
        float rxqidk1      = yr*qidk3 - zr*qidk2;
        float rxqidk2      = zr*qidk1 - xr*qidk3;
        float rxqidk3      = xr*qidk2 - yr*qidk1;
419

420
421
422
        float rxqkdi1      = yr*qkdi3 - zr*qkdi2;
        float rxqkdi2      = zr*qkdi1 - xr*qkdi3;
        float rxqkdi3      = xr*qkdi2 - yr*qkdi1;
423

424
425
426
        float rxqiuk1      = yr*qiuk3 - zr*qiuk2;
        float rxqiuk2      = zr*qiuk1 - xr*qiuk3;
        float rxqiuk3      = xr*qiuk2 - yr*qiuk1;
427

428
429
430
        float rxqkui1      = yr*qkui3 - zr*qkui2;
        float rxqkui2      = zr*qkui1 - xr*qkui3;
        float rxqkui3      = xr*qkui2 - yr*qkui1;
431

432
433
434
        float rxqiukp1     = yr*qiukp3 - zr*qiukp2;
        float rxqiukp2     = zr*qiukp1 - xr*qiukp3;
        float rxqiukp3     = xr*qiukp2 - yr*qiukp1;
435

436
437
438
        float rxqkuip1     = yr*qkuip3 - zr*qkuip2;
        float rxqkuip2     = zr*qkuip1 - xr*qkuip3;
        float rxqkuip3     = xr*qkuip2 - yr*qkuip1;
439
440
441

        // calculate the scalar products for permanent components

442
443
444
445
446
447
448
449
450
        float sc2          = di1*dk1 + di2*dk2 + di3*dk3;
        float sc3          = di1*xr + di2*yr + di3*zr;
        float sc4          = dk1*xr + dk2*yr + dk3*zr;
        float sc5          = qir1*xr + qir2*yr + qir3*zr;
        float sc6          = qkr1*xr + qkr2*yr + qkr3*zr;
        float sc7          = qir1*dk1 + qir2*dk2 + qir3*dk3;
        float sc8          = qkr1*di1 + qkr2*di2 + qkr3*di3;
        float sc9          = qir1*qkr1 + qir2*qkr2 + qir3*qkr3;
        float sc10         = qi1*qk1 + qi2*qk2 + qi3*qk3
451
452
                             + qi4*qk4 + qi5*qk5 + qi6*qk6
                             + qi7*qk7 + qi8*qk8 + qi9*qk9;
453
454
455

        // calculate the scalar products for induced components

456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
        float sci1         = atomI.inducedDipole[0]*dk1 + atomI.inducedDipole[1]*dk2
                             + atomI.inducedDipole[2]*dk3 + di1*atomJ.inducedDipole[0]
                             + di2*atomJ.inducedDipole[1] + di3*atomJ.inducedDipole[2];

        float sci3         = atomI.inducedDipole[0]*xr + atomI.inducedDipole[1]*yr + atomI.inducedDipole[2]*zr;
        float sci4         = atomJ.inducedDipole[0]*xr + atomJ.inducedDipole[1]*yr + atomJ.inducedDipole[2]*zr;
        float sci7         = qir1*atomJ.inducedDipole[0] + qir2*atomJ.inducedDipole[1] + qir3*atomJ.inducedDipole[2];
        float sci8         = qkr1*atomI.inducedDipole[0] + qkr2*atomI.inducedDipole[1] + qkr3*atomI.inducedDipole[2];
        float scip1        = atomI.inducedDipoleP[0]*dk1 + atomI.inducedDipoleP[1]*dk2 + atomI.inducedDipoleP[2]*dk3 + di1*atomJ.inducedDipoleP[0] + di2*atomJ.inducedDipoleP[1] + di3*atomJ.inducedDipoleP[2];
        float scip2        = atomI.inducedDipole[0]*atomJ.inducedDipoleP[0]+atomI.inducedDipole[1]*atomJ.inducedDipoleP[1]
                              + atomI.inducedDipole[2]*atomJ.inducedDipoleP[2]+atomI.inducedDipoleP[0]*atomJ.inducedDipole[0]
                              + atomI.inducedDipoleP[1]*atomJ.inducedDipole[1]+atomI.inducedDipoleP[2]*atomJ.inducedDipole[2];

        float scip3        = atomI.inducedDipoleP[0]*xr + atomI.inducedDipoleP[1]*yr + atomI.inducedDipoleP[2]*zr;
        float scip4        = atomJ.inducedDipoleP[0]*xr + atomJ.inducedDipoleP[1]*yr + atomJ.inducedDipoleP[2]*zr;

        float scip7        = qir1*atomJ.inducedDipoleP[0] + qir2*atomJ.inducedDipoleP[1] + qir3*atomJ.inducedDipoleP[2];
        float scip8        = qkr1*atomI.inducedDipoleP[0] + qkr2*atomI.inducedDipoleP[1] + qkr3*atomI.inducedDipoleP[2];
474
475
476

        // calculate the gl functions for permanent components

477
478
479
480
481
482
483
484
485
        float gl0           = ci*ck;
        float gl1           = ck*sc3 - ci*sc4;
        float gl2           = ci*sc6 + ck*sc5 - sc3*sc4;
        float gl3           = sc3*sc6 - sc4*sc5;
        float gl4           = sc5*sc6;
        float gl5           = -4.0f * sc9;
        float gl6           = sc2;
        float gl7           = 2.0f * (sc7-sc8);
        float gl8           = 2.0f * sc10;
486
487
488

        // calculate the gl functions for induced components

489
490
491
492
493
494
495
496
497
498
        float gli1          = ck*sci3 - ci*sci4;
        float gli2          = -sc3*sci4 - sci3*sc4;
        float gli3          = sci3*sc6 - sci4*sc5;
        float gli6          = sci1;
        float gli7          = 2.0f * (sci7-sci8);
        float glip1         = ck*scip3 - ci*scip4;
        float glip2         = -sc3*scip4 - scip3*sc4;
        float glip3         = scip3*sc6 - scip4*sc5;
        float glip6         = scip1;
        float glip7         = 2.0f * (scip7-scip8);
499
500
501

        // compute the energy contributions for this interaction

502
503
        float e             = bn0*gl0 + bn1*(gl1+gl6) + bn2*(gl2+gl7+gl8) + bn3*(gl3+gl5) + bn4*gl4;
        float ei            = 0.5f * (bn1*(gli1+gli6) + bn2*(gli2+gli7) + bn3*gli3); 
504
505
506

        // get the real energy without any screening function

507
508
509
510
        float erl           = rr1*gl0 + rr3*(gl1+gl6) + rr5*(gl2+gl7+gl8) + rr7*(gl3+gl5) + rr9*gl4;
        float erli          = 0.5f*(rr3*(gli1+gli6)*psc3 + rr5*(gli2+gli7)*psc5 + rr7*gli3*psc7);
        e                   = e - (1.0f-scalingFactors[MScaleIndex])*erl;
        ei                  = ei - erli;
511

Mark Friedrichs's avatar
Mark Friedrichs committed
512
        forceTorqueEnergy[0].w = (e + ei);
513
514
515
516
517
518
519
520

        // increment the total intramolecular energy; assumes;
        // intramolecular distances are less than half of cell;
        // length and less than the ewald cutoff;
/*
        if (molcule(ii) .eq. molcule(kk)) {
           eintra = eintra + mscale(kk)*erl*f;
           eintra = eintra + 0.5f*pscale(kk);
521
522
523
&                        * (rr3*(gli1+gli6)*scale3;
&                              + rr5*(gli2+gli7)*scale5;
&                              + rr7*gli3*scale7);
524
525
526
527
528
        }
*/

        // intermediate variables for permanent force terms

529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
        float gf1 = bn1*gl0 + bn2*(gl1+gl6)
                     + bn3*(gl2+gl7+gl8)
                     + bn4*(gl3+gl5) + bn5*gl4;
        float gf2 = -ck*bn1 + sc4*bn2 - sc6*bn3;
        float gf3 = ci*bn1 + sc3*bn2 + sc5*bn3;
        float gf4 = 2.0f * bn2;
        float gf5 = 2.0f * (-ck*bn2+sc4*bn3-sc6*bn4);
        float gf6 = 2.0f * (-ci*bn2-sc3*bn3-sc5*bn4);
        float gf7 = 4.0f * bn3;
        float gfr1 = rr3*gl0 + rr5*(gl1+gl6)
                      + rr7*(gl2+gl7+gl8)
                      + rr9*(gl3+gl5) + rr11*gl4;
        float gfr2 = -ck*rr3 + sc4*rr5 - sc6*rr7;
        float gfr3 = ci*rr3 + sc3*rr5 + sc5*rr7;
        float gfr4 = 2.0f * rr5;
        float gfr5 = 2.0f * (-ck*rr5+sc4*rr7-sc6*rr9);
        float gfr6 = 2.0f * (-ci*rr5-sc3*rr7-sc5*rr9);
        float gfr7 = 4.0f * rr7;
547
548
549

        // intermediate variables for induced force terms

Peter Eastman's avatar
Peter Eastman committed
550
551
552
553
554
        float gfi1 = 0.5f*(bn2*(gli1+glip1+gli6+glip6)
                      + bn2*scip2
                      + bn3*(gli2+glip2+gli7+glip7)
                      - bn3*(sci3*scip4+scip3*sci4)
                      + bn4*(gli3+glip3));
555
556
557
558
559
        float gfi2 = -ck*bn1 + sc4*bn2 - sc6*bn3;
        float gfi3 = ci*bn1 + sc3*bn2 + sc5*bn3;
        float gfi4 = 2.0f * bn2;
        float gfi5 = bn3 * (sci4+scip4);
        float gfi6 = -bn3 * (sci3+scip3);
Peter Eastman's avatar
Peter Eastman committed
560
        float gfri1 = 0.5f*(rr5*((gli1+gli6)*psc3
561
562
                             + (glip1+glip6)*dsc3
                             + scip2*usc3)
Peter Eastman's avatar
Peter Eastman committed
563
                  + rr7*((gli7+gli2)*psc5
564
565
                             + (glip7+glip2)*dsc5
                      - (sci3*scip4+scip3*sci4)*usc5)
Peter Eastman's avatar
Peter Eastman committed
566
                  + rr9*(gli3*psc7+glip3*dsc7));
567
568
569
        float gfri4 = 2.0f * rr5;
        float gfri5 = rr7 * (sci4*psc7+scip4*dsc7);
        float gfri6 = -rr7 * (sci3*psc7+scip3*dsc7);
570
571
572

        // get the permanent force with screening

573
574
575
576
577
578
579
580
581
        float ftm21 = gf1*xr + gf2*di1 + gf3*dk1
                       + gf4*(qkdi1-qidk1) + gf5*qir1
                       + gf6*qkr1 + gf7*(qiqkr1+qkqir1);
        float ftm22 = gf1*yr + gf2*di2 + gf3*dk2
                       + gf4*(qkdi2-qidk2) + gf5*qir2
                       + gf6*qkr2 + gf7*(qiqkr2+qkqir2);
        float ftm23 = gf1*zr + gf2*di3 + gf3*dk3
                       + gf4*(qkdi3-qidk3) + gf5*qir3
                       + gf6*qkr3 + gf7*(qiqkr3+qkqir3);
582
583
584

        // get the permanent force without screening

585
586
587
588
589
590
591
592
593
        float ftm2r1 = gfr1*xr + gfr2*di1 + gfr3*dk1
                       + gfr4*(qkdi1-qidk1) + gfr5*qir1
                       + gfr6*qkr1 + gfr7*(qiqkr1+qkqir1);
        float ftm2r2 = gfr1*yr + gfr2*di2 + gfr3*dk2
                       + gfr4*(qkdi2-qidk2) + gfr5*qir2
                       + gfr6*qkr2 + gfr7*(qiqkr2+qkqir2);
        float ftm2r3 = gfr1*zr + gfr2*di3 + gfr3*dk3
                       + gfr4*(qkdi3-qidk3) + gfr5*qir3
                       + gfr6*qkr3 + gfr7*(qiqkr3+qkqir3);
594
595
596

        // get the induced force with screening

597
598
599
600
601
602
603
604
605
        float ftm2i1 = gfi1*xr + 0.5f*
              (gfi2*(atomI.inducedDipole[0]+atomI.inducedDipoleP[0])
             + bn2*(sci4*atomI.inducedDipoleP[0]+scip4*atomI.inducedDipole[0])
             + gfi3*(atomJ.inducedDipole[0]+atomJ.inducedDipoleP[0])
             + bn2*(sci3*atomJ.inducedDipoleP[0]+scip3*atomJ.inducedDipole[0])
             + (sci4+scip4)*bn2*di1
             + (sci3+scip3)*bn2*dk1
             + gfi4*(qkui1+qkuip1-qiuk1-qiukp1))
             + gfi5*qir1 + gfi6*qkr1;
Mark Friedrichs's avatar
Mark Friedrichs committed
606

607
608
609
610
611
612
613
614
615
        float ftm2i2 = gfi1*yr + 0.5f*
              (gfi2*(atomI.inducedDipole[1]+atomI.inducedDipoleP[1])
             + bn2*(sci4*atomI.inducedDipoleP[1]+scip4*atomI.inducedDipole[1])
             + gfi3*(atomJ.inducedDipole[1]+atomJ.inducedDipoleP[1])
             + bn2*(sci3*atomJ.inducedDipoleP[1]+scip3*atomJ.inducedDipole[1])
             + (sci4+scip4)*bn2*di2
             + (sci3+scip3)*bn2*dk2
             + gfi4*(qkui2+qkuip2-qiuk2-qiukp2))
             + gfi5*qir2 + gfi6*qkr2;
Mark Friedrichs's avatar
Mark Friedrichs committed
616

617
618
619
620
621
622
623
624
625
        float ftm2i3 = gfi1*zr + 0.5f*
              (gfi2*(atomI.inducedDipole[2]+atomI.inducedDipoleP[2])
             + bn2*(sci4*atomI.inducedDipoleP[2]+scip4*atomI.inducedDipole[2])
             + gfi3*(atomJ.inducedDipole[2]+atomJ.inducedDipoleP[2])
             + bn2*(sci3*atomJ.inducedDipoleP[2]+scip3*atomJ.inducedDipole[2])
             + (sci4+scip4)*bn2*di3
             + (sci3+scip3)*bn2*dk3
             + gfi4*(qkui3+qkuip3-qiuk3-qiukp3))
             + gfi5*qir3 + gfi6*qkr3;
626
627
628

        // get the induced force without screening

629
        float ftm2ri1 = gfri1*xr + 0.5f*
630
            (- rr3*ck*(atomI.inducedDipole[0]*psc3+atomI.inducedDipoleP[0]*dsc3)
631
632
             + rr5*sc4*(atomI.inducedDipole[0]*psc5+atomI.inducedDipoleP[0]*dsc5)
             - rr7*sc6*(atomI.inducedDipole[0]*psc7+atomI.inducedDipoleP[0]*dsc7))
633
             + (rr3*ci*(atomJ.inducedDipole[0]*psc3+atomJ.inducedDipoleP[0]*dsc3)
634
635
636
637
638
639
640
641
642
             + rr5*sc3*(atomJ.inducedDipole[0]*psc5+atomJ.inducedDipoleP[0]*dsc5)
             + rr7*sc5*(atomJ.inducedDipole[0]*psc7+atomJ.inducedDipoleP[0]*dsc7))*0.5f
             + rr5*usc5*(sci4*atomI.inducedDipoleP[0]+scip4*atomI.inducedDipole[0]
             + sci3*atomJ.inducedDipoleP[0]+scip3*atomJ.inducedDipole[0])*0.5f
             + 0.5f*(sci4*psc5+scip4*dsc5)*rr5*di1
             + 0.5f*(sci3*psc5+scip3*dsc5)*rr5*dk1
             + 0.5f*gfri4*((qkui1-qiuk1)*psc5
             + (qkuip1-qiukp1)*dsc5)
             + gfri5*qir1 + gfri6*qkr1;
643

644
        float ftm2ri2 = gfri1*yr + 0.5f*
645
            (- rr3*ck*(atomI.inducedDipole[1]*psc3+atomI.inducedDipoleP[1]*dsc3)
646
647
             + rr5*sc4*(atomI.inducedDipole[1]*psc5+atomI.inducedDipoleP[1]*dsc5)
             - rr7*sc6*(atomI.inducedDipole[1]*psc7+atomI.inducedDipoleP[1]*dsc7))
648
             + (rr3*ci*(atomJ.inducedDipole[1]*psc3+atomJ.inducedDipoleP[1]*dsc3)
649
650
651
652
653
654
655
656
657
             + rr5*sc3*(atomJ.inducedDipole[1]*psc5+atomJ.inducedDipoleP[1]*dsc5)
             + rr7*sc5*(atomJ.inducedDipole[1]*psc7+atomJ.inducedDipoleP[1]*dsc7))*0.5f
             + rr5*usc5*(sci4*atomI.inducedDipoleP[1]+scip4*atomI.inducedDipole[1]
             + sci3*atomJ.inducedDipoleP[1]+scip3*atomJ.inducedDipole[1])*0.5f
             + 0.5f*(sci4*psc5+scip4*dsc5)*rr5*di2
             + 0.5f*(sci3*psc5+scip3*dsc5)*rr5*dk2
             + 0.5f*gfri4*((qkui2-qiuk2)*psc5
             + (qkuip2-qiukp2)*dsc5)
             + gfri5*qir2 + gfri6*qkr2;
658

659
        float ftm2ri3 = gfri1*zr + 0.5f*
660
            (- rr3*ck*(atomI.inducedDipole[2]*psc3+atomI.inducedDipoleP[2]*dsc3)
661
662
             + rr5*sc4*(atomI.inducedDipole[2]*psc5+atomI.inducedDipoleP[2]*dsc5)
             - rr7*sc6*(atomI.inducedDipole[2]*psc7+atomI.inducedDipoleP[2]*dsc7))
663
             + (rr3*ci*(atomJ.inducedDipole[2]*psc3+atomJ.inducedDipoleP[2]*dsc3)
664
665
666
667
668
669
670
671
672
             + rr5*sc3*(atomJ.inducedDipole[2]*psc5+atomJ.inducedDipoleP[2]*dsc5)
             + rr7*sc5*(atomJ.inducedDipole[2]*psc7+atomJ.inducedDipoleP[2]*dsc7))*0.5f
             + rr5*usc5*(sci4*atomI.inducedDipoleP[2]+scip4*atomI.inducedDipole[2]
             + sci3*atomJ.inducedDipoleP[2]+scip3*atomJ.inducedDipole[2])*0.5f
             + 0.5f*(sci4*psc5+scip4*dsc5)*rr5*di3
             + 0.5f*(sci3*psc5+scip3*dsc5)*rr5*dk3
             + 0.5f*gfri4*((qkui3-qiuk3)*psc5
             + (qkuip3-qiukp3)*dsc5)
             + gfri5*qir3 + gfri6*qkr3;
673
674
675

        // account for partially excluded induced interactions

676
677
        float temp3 = 0.5f * rr3 * ((gli1+gli6)*scalingFactors[PScaleIndex]
                                   +(glip1+glip6)*scalingFactors[DScaleIndex]);
Mark Friedrichs's avatar
Mark Friedrichs committed
678

679
680
        float temp5 = 0.5f * rr5 * ((gli2+gli7)*scalingFactors[PScaleIndex]
                                   +(glip2+glip7)*scalingFactors[DScaleIndex]);
Mark Friedrichs's avatar
Mark Friedrichs committed
681

682
683
        float temp7 = 0.5f * rr7 * (gli3*scalingFactors[PScaleIndex]
                                   +glip3*scalingFactors[DScaleIndex]);
Mark Friedrichs's avatar
Mark Friedrichs committed
684

685
686
687
        float fridmp1 = temp3*ddsc31 + temp5*ddsc51 + temp7*ddsc71;
        float fridmp2 = temp3*ddsc32 + temp5*ddsc52 + temp7*ddsc72;
        float fridmp3 = temp3*ddsc33 + temp5*ddsc53 + temp7*ddsc73;
688
689
690

        // find some scaling terms for induced-induced force

Mark Friedrichs's avatar
Mark Friedrichs committed
691
692
        temp3         = 0.5f * rr3 * scalingFactors[UScaleIndex] * scip2;
        temp5         = -0.5f * rr5 * scalingFactors[UScaleIndex] * (sci3*scip4+scip3*sci4);
693
694
695
        float findmp1 = temp3*ddsc31 + temp5*ddsc51;
        float findmp2 = temp3*ddsc32 + temp5*ddsc52;
        float findmp3 = temp3*ddsc33 + temp5*ddsc53;
696
697
698

        // modify the forces for partially excluded interactions

Mark Friedrichs's avatar
Mark Friedrichs committed
699
700
701
        ftm2i1       -= (fridmp1 + findmp1);
        ftm2i2       -= (fridmp2 + findmp2);
        ftm2i3       -= (fridmp3 + findmp3);
702
703
704

        // correction to convert mutual to direct polarization force

705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
        if( cAmoebaSim.polarizationType ){

           float gfd     = 0.5f * (bn2*scip2 - bn3*(scip3*sci4+sci3*scip4));
           ftm2i1       -= gfd*xr + 0.5f*bn2*(sci4*atomI.inducedDipoleP[0]+scip4*atomI.inducedDipole[0]+sci3*atomJ.inducedDipoleP[0]+scip3*atomJ.inducedDipole[0]);
           ftm2i2       -= gfd*yr + 0.5f*bn2*(sci4*atomI.inducedDipoleP[1]+scip4*atomI.inducedDipole[1]+sci3*atomJ.inducedDipoleP[1]+scip3*atomJ.inducedDipole[1]);
           ftm2i3       -= gfd*zr + 0.5f*bn2*(sci4*atomI.inducedDipoleP[2]+scip4*atomI.inducedDipole[2]+sci3*atomJ.inducedDipoleP[2]+scip3*atomJ.inducedDipole[2]);

           float gfdr    = 0.5f * (rr5*scip2*usc3 - rr7*(scip3*sci4 +sci3*scip4)*usc5);
           float fdir1   = gfdr*xr + 0.5f*usc5*rr5*(sci4*atomI.inducedDipoleP[0]+scip4*atomI.inducedDipole[0] + sci3*atomJ.inducedDipoleP[0]+scip3*atomJ.inducedDipole[0]);
           float fdir2   = gfdr*yr + 0.5f*usc5*rr5*(sci4*atomI.inducedDipoleP[1]+scip4*atomI.inducedDipole[1] + sci3*atomJ.inducedDipoleP[1]+scip3*atomJ.inducedDipole[1]);
           float fdir3   = gfdr*zr + 0.5f*usc5*rr5*(sci4*atomI.inducedDipoleP[2]+scip4*atomI.inducedDipole[2] + sci3*atomJ.inducedDipoleP[2]+scip3*atomJ.inducedDipole[2]);

           ftm2i1       += fdir1 + findmp1;
           ftm2i2       += fdir2 + findmp2;
           ftm2i3       += fdir3 + findmp3;
720
721
722
723
        }

        // intermediate variables for induced torque terms

724
725
726
727
728
        float gti2  = 0.5f * bn2 * (sci4+scip4);
        float gti3  = 0.5f * bn2 * (sci3+scip3);
        float gti4  = gfi4;
        float gti5  = gfi5;
        float gti6  = gfi6;
729
730
731
732
733
        float gtri2 = 0.5f * rr5 * (sci4*psc5+scip4*dsc5);
        float gtri3 = 0.5f * rr5 * (sci3*psc5+scip3*dsc5);
        float gtri4 = gfri4;
        float gtri5 = gfri5;
        float gtri6 = gfri6;
734
735
736

        // get the permanent torque with screening

737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
        float ttm21 = -bn1*dixdk1 + gf2*dixr1
            + gf4*(dixqkr1+dkxqir1+rxqidk1-2.0f*qixqk1)
            - gf5*rxqir1 - gf7*(rxqikr1+qkrxqir1);
        float ttm22 = -bn1*dixdk2 + gf2*dixr2
            + gf4*(dixqkr2+dkxqir2+rxqidk2-2.0f*qixqk2)
            - gf5*rxqir2 - gf7*(rxqikr2+qkrxqir2);
        float ttm23 = -bn1*dixdk3 + gf2*dixr3
            + gf4*(dixqkr3+dkxqir3+rxqidk3-2.0f*qixqk3)
            - gf5*rxqir3 - gf7*(rxqikr3+qkrxqir3);
        float ttm31 = bn1*dixdk1 + gf3*dkxr1
            - gf4*(dixqkr1+dkxqir1+rxqkdi1-2.0f*qixqk1)
            - gf6*rxqkr1 - gf7*(rxqkir1-qkrxqir1);
        float ttm32 = bn1*dixdk2 + gf3*dkxr2
            - gf4*(dixqkr2+dkxqir2+rxqkdi2-2.0f*qixqk2)
            - gf6*rxqkr2 - gf7*(rxqkir2-qkrxqir2);
        float ttm33 = bn1*dixdk3 + gf3*dkxr3
            - gf4*(dixqkr3+dkxqir3+rxqkdi3-2.0f*qixqk3)
            - gf6*rxqkr3 - gf7*(rxqkir3-qkrxqir3);
755
756
757

        // get the permanent torque without screening

758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
        float ttm2r1 = -rr3*dixdk1 + gfr2*dixr1-gfr5*rxqir1
            + gfr4*(dixqkr1+dkxqir1+rxqidk1-2.0f*qixqk1)
            - gfr7*(rxqikr1+qkrxqir1);
        float ttm2r2 = -rr3*dixdk2 + gfr2*dixr2-gfr5*rxqir2
            + gfr4*(dixqkr2+dkxqir2+rxqidk2-2.0f*qixqk2)
            - gfr7*(rxqikr2+qkrxqir2);
        float ttm2r3 = -rr3*dixdk3 + gfr2*dixr3-gfr5*rxqir3
            + gfr4*(dixqkr3+dkxqir3+rxqidk3-2.0f*qixqk3)
            - gfr7*(rxqikr3+qkrxqir3);
        float ttm3r1 = rr3*dixdk1 + gfr3*dkxr1 -gfr6*rxqkr1
            - gfr4*(dixqkr1+dkxqir1+rxqkdi1-2.0f*qixqk1)
            - gfr7*(rxqkir1-qkrxqir1);
        float ttm3r2 = rr3*dixdk2 + gfr3*dkxr2 -gfr6*rxqkr2
            - gfr4*(dixqkr2+dkxqir2+rxqkdi2-2.0f*qixqk2)
            - gfr7*(rxqkir2-qkrxqir2);
        float ttm3r3 = rr3*dixdk3 + gfr3*dkxr3 -gfr6*rxqkr3
            - gfr4*(dixqkr3+dkxqir3+rxqkdi3-2.0f*qixqk3)
            - gfr7*(rxqkir3-qkrxqir3);
776
777
778

        // get the induced torque with screening

779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
        float ttm2i1 = -bn1*(dixuk1+dixukp1)*0.5f
            + gti2*dixr1 + gti4*(ukxqir1+rxqiuk1
            + ukxqirp1+rxqiukp1)*0.5f - gti5*rxqir1;
        float ttm2i2 = -bn1*(dixuk2+dixukp2)*0.5f
            + gti2*dixr2 + gti4*(ukxqir2+rxqiuk2
            + ukxqirp2+rxqiukp2)*0.5f - gti5*rxqir2;
        float ttm2i3 = -bn1*(dixuk3+dixukp3)*0.5f
            + gti2*dixr3 + gti4*(ukxqir3+rxqiuk3
            + ukxqirp3+rxqiukp3)*0.5f - gti5*rxqir3;
        float ttm3i1 = -bn1*(dkxui1+dkxuip1)*0.5f
            + gti3*dkxr1 - gti4*(uixqkr1+rxqkui1
            + uixqkrp1+rxqkuip1)*0.5f - gti6*rxqkr1;
        float ttm3i2 = -bn1*(dkxui2+dkxuip2)*0.5f
            + gti3*dkxr2 - gti4*(uixqkr2+rxqkui2
            + uixqkrp2+rxqkuip2)*0.5f - gti6*rxqkr2;
        float ttm3i3 = -bn1*(dkxui3+dkxuip3)*0.5f
            + gti3*dkxr3 - gti4*(uixqkr3+rxqkui3
            + uixqkrp3+rxqkuip3)*0.5f - gti6*rxqkr3;
797
798
799

        // get the induced torque without screening

800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
        float ttm2ri1 = -rr3*(dixuk1*psc3+dixukp1*dsc3)*0.5f
            + gtri2*dixr1 + gtri4*((ukxqir1+rxqiuk1)*psc5
            +(ukxqirp1+rxqiukp1)*dsc5)*0.5f - gtri5*rxqir1;
        float ttm2ri2 = -rr3*(dixuk2*psc3+dixukp2*dsc3)*0.5f
            + gtri2*dixr2 + gtri4*((ukxqir2+rxqiuk2)*psc5
            +(ukxqirp2+rxqiukp2)*dsc5)*0.5f - gtri5*rxqir2;
        float ttm2ri3 = -rr3*(dixuk3*psc3+dixukp3*dsc3)*0.5f
            + gtri2*dixr3 + gtri4*((ukxqir3+rxqiuk3)*psc5
            +(ukxqirp3+rxqiukp3)*dsc5)*0.5f - gtri5*rxqir3;
        float ttm3ri1 = -rr3*(dkxui1*psc3+dkxuip1*dsc3)*0.5f
            + gtri3*dkxr1 - gtri4*((uixqkr1+rxqkui1)*psc5
            +(uixqkrp1+rxqkuip1)*dsc5)*0.5f - gtri6*rxqkr1;
        float ttm3ri2 = -rr3*(dkxui2*psc3+dkxuip2*dsc3)*0.5f
            + gtri3*dkxr2 - gtri4*((uixqkr2+rxqkui2)*psc5
            +(uixqkrp2+rxqkuip2)*dsc5)*0.5f - gtri6*rxqkr2;
        float ttm3ri3 = -rr3*(dkxui3*psc3+dkxuip3*dsc3)*0.5f
            + gtri3*dkxr3 - gtri4*((uixqkr3+rxqkui3)*psc5
            +(uixqkrp3+rxqkuip3)*dsc5)*0.5f - gtri6*rxqkr3;
818
819
820

        // handle the case where scaling is used

821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
        ftm21  = (ftm21-(1.0f-scalingFactors[MScaleIndex])*ftm2r1);
        ftm2i1 = (ftm2i1-ftm2ri1);
        ttm21  = (ttm21-(1.0f-scalingFactors[MScaleIndex])*ttm2r1);
        ttm2i1 = (ttm2i1-ttm2ri1);
        ttm31  = (ttm31-(1.0f-scalingFactors[MScaleIndex])*ttm3r1);
        ttm3i1 = (ttm3i1-ttm3ri1);

        ftm22  = (ftm22-(1.0f-scalingFactors[MScaleIndex])*ftm2r2);
        ftm2i2 = (ftm2i2-ftm2ri2);
        ttm22  = (ttm22-(1.0f-scalingFactors[MScaleIndex])*ttm2r2);
        ttm2i2 = (ttm2i2-ttm2ri2);
        ttm32  = (ttm32-(1.0f-scalingFactors[MScaleIndex])*ttm3r2);
        ttm3i2 = (ttm3i2-ttm3ri2);

        ftm23  = (ftm23-(1.0f-scalingFactors[MScaleIndex])*ftm2r3);
        ftm2i3 = (ftm2i3-ftm2ri3);
        ttm23  = (ttm23-(1.0f-scalingFactors[MScaleIndex])*ttm2r3);
        ttm2i3 = (ttm2i3-ttm2ri3);
        ttm33  = (ttm33-(1.0f-scalingFactors[MScaleIndex])*ttm3r3);
        ttm3i3 = (ttm3i3-ttm3ri3);
841
842

        // increment gradient due to force and torque on first site;
Mark Friedrichs's avatar
Mark Friedrichs committed
843

Mark Friedrichs's avatar
Mark Friedrichs committed
844
845
846
        forceTorqueEnergy[0].x       = (ftm21 + ftm2i1);
        forceTorqueEnergy[0].y       = (ftm22 + ftm2i2);
        forceTorqueEnergy[0].z       = (ftm23 + ftm2i3);
Peter Eastman's avatar
Peter Eastman committed
847

Mark Friedrichs's avatar
Mark Friedrichs committed
848
849
850
        forceTorqueEnergy[1].x       =  (ttm21 + ttm2i1);
        forceTorqueEnergy[1].y       =  (ttm22 + ttm2i2);
        forceTorqueEnergy[1].z       =  (ttm23 + ttm2i3);
Peter Eastman's avatar
Peter Eastman committed
851

Mark Friedrichs's avatar
Mark Friedrichs committed
852
853
854
        forceTorqueEnergy[2].x       =  (ttm31 + ttm3i1);
        forceTorqueEnergy[2].y       =  (ttm32 + ttm3i2);
        forceTorqueEnergy[2].z       =  (ttm33 + ttm3i3);
Mark Friedrichs's avatar
Mark Friedrichs committed
855
856
857
858

#ifdef AMOEBA_DEBUG
    int debugIndex               = 0;
    float idTracker              = 1.0f;
859
/*
Mark Friedrichs's avatar
Mark Friedrichs committed
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
    debugArray[debugIndex].x = atomI.labFrameDipole[0];
    debugArray[debugIndex].y = atomI.labFrameDipole[1];
    debugArray[debugIndex].z = atomI.labFrameDipole[2];
    debugArray[debugIndex].w = r2;

    debugIndex++;
    idTracker               += 1.0;
    debugArray[debugIndex].x = atomJ.labFrameDipole[0];
    debugArray[debugIndex].y = atomJ.labFrameDipole[1];
    debugArray[debugIndex].z = atomJ.labFrameDipole[2];
    debugArray[debugIndex].w = cSim.alphaEwald;

    debugIndex++;
    idTracker               += 1.0;
    debugArray[debugIndex].x = atomI.inducedDipole[0];
    debugArray[debugIndex].y = atomI.inducedDipole[1];
    debugArray[debugIndex].z = atomI.inducedDipole[2];
    debugArray[debugIndex].w = idTracker;

    debugIndex++;
    idTracker               += 1.0;
    debugArray[debugIndex].x = atomJ.inducedDipole[0];
    debugArray[debugIndex].y = atomJ.inducedDipole[1];
    debugArray[debugIndex].z = atomJ.inducedDipole[2];
    debugArray[debugIndex].w = idTracker;

    debugIndex++;
    idTracker               += 1.0;
    debugArray[debugIndex].x = atomI.inducedDipoleP[0];
    debugArray[debugIndex].y = atomI.inducedDipoleP[1];
    debugArray[debugIndex].z = atomI.inducedDipoleP[2];
    debugArray[debugIndex].w = idTracker;

    debugIndex++;
    idTracker               += 1.0;
    debugArray[debugIndex].x = atomJ.inducedDipoleP[0];
    debugArray[debugIndex].y = atomJ.inducedDipoleP[1];
    debugArray[debugIndex].z = atomJ.inducedDipoleP[2];
    debugArray[debugIndex].w = idTracker;

    debugIndex++;
    idTracker               += 1.0;
902
903
904
    debugArray[debugIndex].x = conversionFactor*ftm21;
    debugArray[debugIndex].y = conversionFactor*ftm22;
    debugArray[debugIndex].z = conversionFactor*ftm23;
Mark Friedrichs's avatar
Mark Friedrichs committed
905
906
    debugArray[debugIndex].w = idTracker;
    debugIndex++;
Mark Friedrichs's avatar
Mark Friedrichs committed
907

Mark Friedrichs's avatar
Mark Friedrichs committed
908
    idTracker               += 1.0;
909
910
911
912
    debugArray[debugIndex].x = e;
    debugArray[debugIndex].y = ei;
    debugArray[debugIndex].z = erl;
    debugArray[debugIndex].w = erli;
Mark Friedrichs's avatar
Mark Friedrichs committed
913
    debugIndex++;
Mark Friedrichs's avatar
Mark Friedrichs committed
914
915
916
*/

    idTracker               += 100.0;
917
918
    debugArray[debugIndex].x = r2;
    debugArray[debugIndex].y = cSim.alphaEwald;
Mark Friedrichs's avatar
Mark Friedrichs committed
919
920
921
    debugArray[debugIndex].z = conversionFactor;
    debugArray[debugIndex].w = idTracker;
    debugIndex++;
Mark Friedrichs's avatar
Mark Friedrichs committed
922

Mark Friedrichs's avatar
Mark Friedrichs committed
923
924
925
926
927
    idTracker               += 100.0;
    debugArray[debugIndex].x = conversionFactor*ftm21;
    debugArray[debugIndex].y = conversionFactor*ftm22;
    debugArray[debugIndex].z = conversionFactor*ftm23;
    debugArray[debugIndex].w = idTracker;
Mark Friedrichs's avatar
Mark Friedrichs committed
928
929
    debugIndex++;

Mark Friedrichs's avatar
Mark Friedrichs committed
930
931
932
933
934
    idTracker               += 100.0;
    debugArray[debugIndex].x = conversionFactor*ftm2i1;
    debugArray[debugIndex].y = conversionFactor*ftm2i2;
    debugArray[debugIndex].z = conversionFactor*ftm2i3;
    debugArray[debugIndex].w = idTracker;
Mark Friedrichs's avatar
Mark Friedrichs committed
935
    debugIndex++;
Mark Friedrichs's avatar
Mark Friedrichs committed
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
/*
    idTracker               += 100.0;
    debugArray[debugIndex].x = fridmp1;
    debugArray[debugIndex].y = fridmp2;
    debugArray[debugIndex].z = fridmp3;
    debugArray[debugIndex].w = idTracker;
    debugIndex++;

    idTracker               += 100.0;
    debugArray[debugIndex].x = findmp1;
    debugArray[debugIndex].y = findmp2;
    debugArray[debugIndex].z = findmp3;
    debugArray[debugIndex].w = idTracker;
    debugIndex++;
*/

    idTracker               += 100.0;
953
954
955
    debugArray[debugIndex].x = conversionFactor*ttm21;
    debugArray[debugIndex].y = conversionFactor*ttm22;
    debugArray[debugIndex].z = conversionFactor*ttm23;
Mark Friedrichs's avatar
Mark Friedrichs committed
956
957
    debugArray[debugIndex].w = idTracker;
    debugIndex++;
Mark Friedrichs's avatar
Mark Friedrichs committed
958
959

    idTracker               += 100.0;
960
961
962
    debugArray[debugIndex].x = conversionFactor*ttm2i1;
    debugArray[debugIndex].y = conversionFactor*ttm2i2;
    debugArray[debugIndex].z = conversionFactor*ttm2i3;
Mark Friedrichs's avatar
Mark Friedrichs committed
963
    debugArray[debugIndex].w = idTracker;
Mark Friedrichs's avatar
Mark Friedrichs committed
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
    debugIndex++;

    idTracker               += 100.0;
    debugArray[debugIndex].x = conversionFactor*ttm31;
    debugArray[debugIndex].y = conversionFactor*ttm32;
    debugArray[debugIndex].z = conversionFactor*ttm33;
    debugArray[debugIndex].w = idTracker;
    debugIndex++;

    idTracker               += 100.0;
    debugArray[debugIndex].x = conversionFactor*ttm3i1;
    debugArray[debugIndex].y = conversionFactor*ttm3i2;
    debugArray[debugIndex].z = conversionFactor*ttm3i3;
    debugArray[debugIndex].w = idTracker;
    debugIndex++;

Mark Friedrichs's avatar
Mark Friedrichs committed
980
#endif
981

Mark Friedrichs's avatar
Mark Friedrichs committed
982
983
    } else {

Mark Friedrichs's avatar
Mark Friedrichs committed
984
985
986
        forceTorqueEnergy[0].x  = 0.0f;
        forceTorqueEnergy[0].y  = 0.0f;
        forceTorqueEnergy[0].z  = 0.0f;
Peter Eastman's avatar
Peter Eastman committed
987

Mark Friedrichs's avatar
Mark Friedrichs committed
988
989
990
        forceTorqueEnergy[1].x  = 0.0f;
        forceTorqueEnergy[1].y  = 0.0f;
        forceTorqueEnergy[1].z  = 0.0f;
Peter Eastman's avatar
Peter Eastman committed
991

Mark Friedrichs's avatar
Mark Friedrichs committed
992
993
994
        forceTorqueEnergy[2].x  = 0.0f;
        forceTorqueEnergy[2].y  = 0.0f;
        forceTorqueEnergy[2].z  = 0.0f;
Mark Friedrichs's avatar
Mark Friedrichs committed
995

Mark Friedrichs's avatar
Mark Friedrichs committed
996
        forceTorqueEnergy[0].w  = 0.0f;
997

Mark Friedrichs's avatar
Mark Friedrichs committed
998
#ifdef AMOEBA_DEBUG
Mark Friedrichs's avatar
Mark Friedrichs committed
999
for( int ii = 0; ii < 12; ii++ ){
Mark Friedrichs's avatar
Mark Friedrichs committed
1000
1001
1002
    debugArray[ii].x = 0.0f;
    debugArray[ii].y = 0.0f;
    debugArray[ii].z = 0.0f;
Mark Friedrichs's avatar
Mark Friedrichs committed
1003
    debugArray[ii].w = (float) (-ii);
Mark Friedrichs's avatar
Mark Friedrichs committed
1004
1005
}
#endif
Peter Eastman's avatar
Peter Eastman committed
1006

1007
1008
1009
1010
1011
1012
    }

    return;

}

Mark Friedrichs's avatar
Mark Friedrichs committed
1013
__device__ void loadPmeDirectElectrostaticShared( struct PmeDirectElectrostaticParticle* sA, unsigned int atomI )
1014
1015
{
    // coordinates & charge
1016
1017
1018
1019
1020
    float4 posq                  = cSim.pPosq[atomI];
    sA->x                        = posq.x;
    sA->y                        = posq.y;
    sA->z                        = posq.z;
    sA->q                        = posq.w;
1021
1022
1023

    // lab dipole

Mark Friedrichs's avatar
Mark Friedrichs committed
1024
1025
1026
1027
    sA->labFrameDipole[0]        = cAmoebaSim.pLabFrameDipole[atomI*3];
    sA->labFrameDipole[1]        = cAmoebaSim.pLabFrameDipole[atomI*3+1];
    sA->labFrameDipole[2]        = cAmoebaSim.pLabFrameDipole[atomI*3+2];

1028
1029
1030

    // lab quadrupole

Mark Friedrichs's avatar
Mark Friedrichs committed
1031
1032
1033
1034
1035
1036
1037
1038
1039
    sA->labFrameQuadrupole[0]    = cAmoebaSim.pLabFrameQuadrupole[atomI*9];
    sA->labFrameQuadrupole[1]    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+1];
    sA->labFrameQuadrupole[2]    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+2];
    sA->labFrameQuadrupole[3]    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+3];
    sA->labFrameQuadrupole[4]    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+4];
    sA->labFrameQuadrupole[5]    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+5];
    sA->labFrameQuadrupole[6]    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+6];
    sA->labFrameQuadrupole[7]    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+7];
    sA->labFrameQuadrupole[8]    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+8];
1040
1041
1042

    // induced dipole

Mark Friedrichs's avatar
Mark Friedrichs committed
1043
1044
1045
    sA->inducedDipole[0]         =  cAmoebaSim.pInducedDipole[atomI*3];
    sA->inducedDipole[1]         =  cAmoebaSim.pInducedDipole[atomI*3+1];
    sA->inducedDipole[2]         =  cAmoebaSim.pInducedDipole[atomI*3+2];
1046
1047
1048

    // induced dipole polar

Mark Friedrichs's avatar
Mark Friedrichs committed
1049
1050
1051
    sA->inducedDipoleP[0]        =  cAmoebaSim.pInducedDipolePolar[atomI*3];
    sA->inducedDipoleP[1]        =  cAmoebaSim.pInducedDipolePolar[atomI*3+1];
    sA->inducedDipoleP[2]        =  cAmoebaSim.pInducedDipolePolar[atomI*3+2];
1052

1053
1054
1055
    float2 dampingFactorAndThole = cAmoebaSim.pDampingFactorAndThole[atomI];
    sA->damp                     = dampingFactorAndThole.x;
    sA->thole                    = dampingFactorAndThole.y;
1056
1057
1058
1059
1060
1061

}

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

#undef USE_OUTPUT_BUFFER_PER_WARP
Mark Friedrichs's avatar
Mark Friedrichs committed
1062
#define METHOD_NAME(a, b) a##Cutoff##b
Mark Friedrichs's avatar
Mark Friedrichs committed
1063
#include "kCalculateAmoebaCudaPmeDirectElectrostatic.h"
1064
1065
#define USE_OUTPUT_BUFFER_PER_WARP
#undef METHOD_NAME
Mark Friedrichs's avatar
Mark Friedrichs committed
1066
#define METHOD_NAME(a, b) a##CutoffByWarp##b
Mark Friedrichs's avatar
Mark Friedrichs committed
1067
#include "kCalculateAmoebaCudaPmeDirectElectrostatic.h"
1068

1069
// reduce psWorkArray_3_1 -> torque
1070

1071
static void kReduceTorque(amoebaGpuContext amoebaGpu )
1072
{
1073
1074
    gpuContext gpu = amoebaGpu->gpuContext;
    kReduceFields_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.bsf_reduce_threads_per_block>>>(
1075
                           gpu->sim.paddedNumberOfAtoms*3, gpu->sim.outputBuffers,
1076
                           amoebaGpu->psWorkArray_3_1->_pDevData, amoebaGpu->psTorque->_pDevData, 1 );
Mark Friedrichs's avatar
Mark Friedrichs committed
1077
    LAUNCHERROR("kReducePmeDirectElectrostaticTorque");
1078
1079
1080
1081
}

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

Mark Friedrichs's avatar
Mark Friedrichs committed
1082
   Compute Amoeba dirrect space portion of electrostatic force & torque
1083
1084
1085
1086
1087

   @param amoebaGpu        amoebaGpu context

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

Mark Friedrichs's avatar
Mark Friedrichs committed
1088
void cudaComputeAmoebaPmeDirectElectrostatic( amoebaGpuContext amoebaGpu )
1089
{
Peter Eastman's avatar
Peter Eastman committed
1090

1091
#ifdef AMOEBA_DEBUG
Mark Friedrichs's avatar
Mark Friedrichs committed
1092
    static const char* methodName = "cudaComputeAmoebaPmeDirectElectrostatic";
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
    static int timestep = 0;
    std::vector<int> fileId;
    timestep++;
    fileId.resize( 2 );
    fileId[0] = timestep;
    fileId[1] = 1;
#endif

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

    gpuContext gpu = amoebaGpu->gpuContext;

    // apparently debug array can take up nontrivial no. registers

#ifdef AMOEBA_DEBUG
    if( amoebaGpu->log ){
Mark Friedrichs's avatar
Mark Friedrichs committed
1109
1110
      (void) fprintf( amoebaGpu->log, "%s %d maxCovalentDegreeSz=%d ZZZ\n",
                      methodName, gpu->natoms, amoebaGpu->maxCovalentDegreeSz );
Peter Eastman's avatar
Peter Eastman committed
1111
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
1112
    int paddedNumberOfAtoms                   = amoebaGpu->gpuContext->sim.paddedNumberOfAtoms;
Mark Friedrichs's avatar
Mark Friedrichs committed
1113
1114
1115
    int maxOffset                             = 20;
    CUDAStream<float4>* debugArray            = new CUDAStream<float4>(maxOffset*paddedNumberOfAtoms, 1, "DebugArray");
    memset( debugArray->_pSysData,      0, sizeof( float )*4*maxOffset*paddedNumberOfAtoms);
1116
    debugArray->Upload();
Mark Friedrichs's avatar
Mark Friedrichs committed
1117
    unsigned int targetAtom                   = 49;
Mark Friedrichs's avatar
Mark Friedrichs committed
1118
1119
#endif

1120
1121
    // on first pass, set threads/block

1122
    static unsigned int threadsPerBlock = 0;
1123
    if( threadsPerBlock == 0 ){
Mark Friedrichs's avatar
Mark Friedrichs committed
1124
1125
1126
1127
1128
1129
1130
        unsigned int maxThreads;
        if (gpu->sm_version >= SM_20)
            maxThreads = 384;
        else if (gpu->sm_version >= SM_12)
            maxThreads = 128;
        else
            maxThreads = 64;
Mark Friedrichs's avatar
Mark Friedrichs committed
1131
        threadsPerBlock = std::min(getThreadsPerBlock(amoebaGpu, sizeof(PmeDirectElectrostaticParticle)), maxThreads);
1132
1133
    }

1134
    kClearFields_3( amoebaGpu, 1 );
1135

Mark Friedrichs's avatar
Mark Friedrichs committed
1136
#ifdef AMOEBA_DEBUG
1137
    if( amoebaGpu->log ){
Mark Friedrichs's avatar
Mark Friedrichs committed
1138
        (void) fprintf( amoebaGpu->log, "kCalculateAmoebaPmeDirectElectrostaticCutoffForces:  numBlocks=%u numThreads=%u bufferPerWarp=%u atm=%lu shrd=%lu ixnCt=%lu workUnits=%u gpu->nonbond_threads_per_block=%u\n",
1139
1140
1141
1142
                        gpu->sim.nonbond_blocks, threadsPerBlock, gpu->bOutputBufferPerWarp,
                        sizeof(PmeDirectElectrostaticParticle), (sizeof(PmeDirectElectrostaticParticle))*threadsPerBlock,
                        (*gpu->psInteractionCount)[0], gpu->sim.workUnits, gpu->sim.nonbond_threads_per_block );
        (void) fflush( amoebaGpu->log );
1143
    }   
Mark Friedrichs's avatar
Mark Friedrichs committed
1144
1145
#endif

1146
1147
    if (gpu->bOutputBufferPerWarp){

1148
        kCalculateAmoebaPmeDirectElectrostaticCutoffByWarpForces_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock, sizeof(PmeDirectElectrostaticParticle)*threadsPerBlock>>>(
Mark Friedrichs's avatar
Mark Friedrichs committed
1149
                                                                          gpu->sim.pInteractingWorkUnit,
1150
#ifdef AMOEBA_DEBUG
1151
                                                                          amoebaGpu->psWorkArray_3_1->_pDevData,
Mark Friedrichs's avatar
Mark Friedrichs committed
1152
                                                                          debugArray->_pDevData, targetAtom );
1153
#else
1154
                                                                          amoebaGpu->psWorkArray_3_1->_pDevData );
1155
1156
1157
1158
#endif

    } else {

Mark Friedrichs's avatar
Mark Friedrichs committed
1159
1160
1161
1162
1163
1164
/*
        if (gpu->sm_version >= SM_20)
            cudaFuncSetCacheConfig(kCalculateAmoebaPmeDirectElectrostaticCutoffForces_kernel, cudaFuncCachePreferL1 );
            //cudaFuncSetCacheConfig(kCalculateAmoebaPmeDirectElectrostaticCutoffForces_kernel, cudaFuncCachePreferShared );
*/

1165
        kCalculateAmoebaPmeDirectElectrostaticCutoffForces_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock, sizeof(PmeDirectElectrostaticParticle)*threadsPerBlock>>>(
Mark Friedrichs's avatar
Mark Friedrichs committed
1166
                                                                    gpu->sim.pInteractingWorkUnit,
1167
#ifdef AMOEBA_DEBUG
1168
                                                                    amoebaGpu->psWorkArray_3_1->_pDevData,
Mark Friedrichs's avatar
Mark Friedrichs committed
1169
                                                                    debugArray->_pDevData, targetAtom );
1170
#else
1171
                                                                    amoebaGpu->psWorkArray_3_1->_pDevData );
1172
1173
#endif
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
1174
    LAUNCHERROR("kCalculateAmoebaPmeDirectElectrostaticCutoffForces");
1175

1176
    kReduceTorque( amoebaGpu );
Mark Friedrichs's avatar
Mark Friedrichs committed
1177

Mark Friedrichs's avatar
Mark Friedrichs committed
1178
}
1179

Mark Friedrichs's avatar
Mark Friedrichs committed
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
/**---------------------------------------------------------------------------------------

   Compute Amoeba electrostatic force & torque using PME

   @param amoebaGpu        amoebaGpu context

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

void cudaComputeAmoebaPmeElectrostatic( amoebaGpuContext amoebaGpu )
{
1190

Mark Friedrichs's avatar
Mark Friedrichs committed
1191
    cudaComputeAmoebaPmeDirectElectrostatic( amoebaGpu );
1192
    kCalculateAmoebaPMEInducedDipoleForces( amoebaGpu );
1193
    cudaComputeAmoebaMapTorqueAndAddToForce( amoebaGpu, amoebaGpu->psTorque );
Peter Eastman's avatar
Peter Eastman committed
1194
}
Mark Friedrichs's avatar
Mark Friedrichs committed
1195