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

#include <stdio.h>
#include <cuda.h>
#include <vector_functions.h>
#include <cstdlib>
#include <string>
#include <iostream>
#include <fstream>
using namespace std;

#include "gputypes.h"

extern __shared__ Vectors sV[];
static __constant__ cudaGmxSimulation cSim;

41
42
/* Cuda compiler on Windows does not recognized "static const float" values */
#define LOCAL_HACK_PI 3.1415926535897932384626433832795
Peter Eastman's avatar
Peter Eastman committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

#define DOT3(v1, v2) (v1.x * v2.x + v1.y * v2.y + v1.z * v2.z)

#define GETNORMEDDOTPRODUCT(v1, v2, dp) \
{ \
    dp          = DOT3(v1, v2); \
    float norm1 = DOT3(v1, v1); \
    float norm2 = DOT3(v2, v2); \
    dp /= sqrt(norm1 * norm2); \
    dp = min(dp, 1.0f); \
    dp = max(dp, -1.0f); \
}

#define CROSS_PRODUCT(v1, v2, c) \
    c.x = v1.y * v2.z - v1.z * v2.y; \
    c.y = v1.z * v2.x - v1.x * v2.z; \
    c.z = v1.x * v2.y - v1.y * v2.x;

#define GETPREFACTORSGIVENANGLECOSINE(cosine, param, dEdR) \
{ \
   float angle          = acos(cosine); \
64
   float deltaIdeal     = angle - (param.x * (LOCAL_HACK_PI / 180.0f)); \
Peter Eastman's avatar
Peter Eastman committed
65
66
67
   dEdR                 = param.y * deltaIdeal; \
}

68
69
70
#define GETENERGYGIVENANGLECOSINE(cosine, param, dEdR) \
{ \
   float angle          = acos(cosine); \
71
   float deltaIdeal     = angle - (param.x * (LOCAL_HACK_PI / 180.0f)); \
72
73
74
   dEdR                 = param.y * deltaIdeal * deltaIdeal; \
}

Peter Eastman's avatar
Peter Eastman committed
75
76
77
78
#define GETANGLEBETWEENTWOVECTORS(v1, v2, angle) \
{ \
    float dp; \
    GETNORMEDDOTPRODUCT(v1, v2, dp); \
79
80
81
82
83
84
85
86
87
88
89
    if (dp > 0.99f || dp < -0.99f) { \
        float4 cross; \
        CROSS_PRODUCT(v1, v2, cross); \
        float scale = DOT3(v1, v1)*DOT3(v2, v2); \
        angle = asin(sqrt(DOT3(cross, cross)/scale)); \
        if (dp < 0.0f) \
            angle = LOCAL_HACK_PI-angle; \
    } \
    else { \
        angle = acos(dp); \
    } \
Peter Eastman's avatar
Peter Eastman committed
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
}

#define GETANGLECOSINEBETWEENTWOVECTORS(v1, v2, angle, cosine) \
{ \
    GETNORMEDDOTPRODUCT(v1, v2, cosine); \
    angle = acos(cosine); \
}

#define GETDIHEDRALANGLEBETWEENTHREEVECTORS(vector1, vector2, vector3, signVector, cp0, cp1, angle) \
{ \
    CROSS_PRODUCT(vector1, vector2, cp0); \
    CROSS_PRODUCT(vector2, vector3, cp1); \
    GETANGLEBETWEENTWOVECTORS(cp0, cp1, angle); \
    float dp = DOT3(signVector, cp1); \
    angle = (dp >= 0) ? angle : -angle; \
}                                                          

#define GETDIHEDRALANGLECOSINEBETWEENTHREEVECTORS(vector1, vector2, vector3, signVector, cp0, cp1, angle, cosine) \
{ \
    CROSS_PRODUCT(vector1, vector2, cp0); \
    CROSS_PRODUCT(vector2, vector3, cp1); \
    GETANGLECOSINEBETWEENTWOVECTORS(cp0, cp1, angle, cosine); \
    float dp = DOT3(signVector, cp1); \
    angle = (dp >= 0) ? angle : -angle; \
}

void SetCalculateLocalForcesSim(gpuContext gpu)
{
    cudaError_t status;
    status = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));     
    RTERROR(status, "cudaMemcpyToSymbol: SetSim copy to cSim failed");
}

void GetCalculateLocalForcesSim(gpuContext gpu)
{
    cudaError_t status;
    status = cudaMemcpyFromSymbol(&gpu->sim, cSim, sizeof(cudaGmxSimulation));     
    RTERROR(status, "cudaMemcpyFromSymbol: SetSim copy from cSim failed");
}
    

131
__global__ void kCalculateLocalForces_kernel()
Peter Eastman's avatar
Peter Eastman committed
132
133
134
135
{
    unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;
    Vectors* A = &sV[threadIdx.x];

136
137
    float energy = 0.0f;

Peter Eastman's avatar
Peter Eastman committed
138
139
140
141
142
143
144
145
146
147
148
149
150
151
    while (pos < cSim.bond_offset)
    {
        if (pos < cSim.bonds)
        {
            int4   atom         = cSim.pBondID[pos];
            float4 atomA        = cSim.pPosq[atom.x];
            float4 atomB        = cSim.pPosq[atom.y];
            float2 bond         = cSim.pBondParameter[pos];
            float dx            = atomB.x - atomA.x;
            float dy            = atomB.y - atomA.y;
            float dz            = atomB.z - atomA.z;
            float r2            = dx * dx + dy * dy + dz * dz;
            float r             = sqrt(r2);
            float deltaIdeal    = r - bond.x;
152
/* E */     energy             += 0.5f * bond.y * deltaIdeal * deltaIdeal;
Peter Eastman's avatar
Peter Eastman committed
153
154
155
156
157
158
159
160
            float dEdR          = bond.y * deltaIdeal;
            dEdR                = (r > 0.0f) ? (dEdR / r) : 0.0f;
//            printf("D: %11.4f %11.4f %11.4f %11.4f %11.4f %11.4f\n", dx, dy, dz, r, deltaIdeal, dEdR);
            dx                 *= dEdR;
            dy                 *= dEdR;
            dz                 *= dEdR;
            unsigned int offsetA                = atom.x + atom.z * cSim.stride;
            unsigned int offsetB                = atom.y + atom.w * cSim.stride;
161
162
            float4 forceA                       = cSim.pForce4[offsetA];
            float4 forceB                       = cSim.pForce4[offsetB];
Peter Eastman's avatar
Peter Eastman committed
163
164
165
166
167
168
169
170
171
172
173
            forceA.x                           += dx;
            forceA.y                           += dy;
            forceA.z                           += dz;
            forceB.x                           -= dx;
            forceB.y                           -= dy;
            forceB.z                           -= dz;
            cSim.pForce4[offsetA]               = forceA;
            cSim.pForce4[offsetB]               = forceB;    
        }
        pos += blockDim.x * gridDim.x;
    }
174

Peter Eastman's avatar
Peter Eastman committed
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
    while (pos < cSim.bond_angle_offset)
    {
        unsigned int pos1   = pos - cSim.bond_offset;
        if (pos1 < cSim.bond_angles)
        {
            int4   atom1            = cSim.pBondAngleID1[pos1];  
            float2 bond_angle       = cSim.pBondAngleParameter[pos1];
            float4 a1               = cSim.pPosq[atom1.x];
            float4 a2               = cSim.pPosq[atom1.y];
            float4 a3               = cSim.pPosq[atom1.z];
            A->v0.x                 = a2.x - a1.x;
            A->v0.y                 = a2.y - a1.y;
            A->v0.z                 = a2.z - a1.z;
            A->v1.x                 = a2.x - a3.x;
            A->v1.y                 = a2.y - a3.y;
            A->v1.z                 = a2.z - a3.z;
            float3 cp;
            CROSS_PRODUCT(A->v0, A->v1, cp);
            float rp                = DOT3(cp, cp); //cx * cx + cy * cy + cz * cz;
            rp                      = max(sqrt(rp), 1.0e-06f);
            float r21               = DOT3(A->v0, A->v0); // dx1 * dx1 + dy1 * dy1 + dz1 * dz1;
            float r23               = DOT3(A->v1, A->v1); // dx2 * dx2 + dy2 * dy2 + dz2 * dz2;
            float dot               = DOT3(A->v0, A->v1); // dx1 * dx2 + dy1 * dy2 + dz1 * dz2;
            float cosine            = dot / sqrt(r21 * r23);
199
200
201

            float angle_energy;
/* E */     GETENERGYGIVENANGLECOSINE(cosine, bond_angle, angle_energy);
202
            energy                 += 0.5f*angle_energy;
203

Peter Eastman's avatar
Peter Eastman committed
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
            float dEdR;
            GETPREFACTORSGIVENANGLECOSINE(cosine, bond_angle, dEdR);
            //printf("%11.4f %11.4f\n", cosine, dEdR);
            float termA             =  dEdR / (r21 * rp);
            float termC             = -dEdR / (r23 * rp);
            float3 c21;
            float3 c23;
            CROSS_PRODUCT(A->v0, cp, c21);
            CROSS_PRODUCT(A->v1, cp, c23);
            c21.x                  *= termA;
            c21.y                  *= termA;
            c21.z                  *= termA;
            c23.x                  *= termC;
            c23.y                  *= termC;
            c23.z                  *= termC;
            int2 atom2              = cSim.pBondAngleID2[pos1];
            unsigned int offset     = atom1.x + atom1.w * cSim.stride;
221
            float4 force            = cSim.pForce4[offset]; 
Peter Eastman's avatar
Peter Eastman committed
222
223
224
225
226
            force.x                += c21.x;
            force.y                += c21.y;
            force.z                += c21.z;
            cSim.pForce4[offset]    = force;
            offset                  = atom1.y + atom2.x * cSim.stride;
227
            force                   = cSim.pForce4[offset];
Peter Eastman's avatar
Peter Eastman committed
228
229
230
231
232
            force.x                -= (c21.x + c23.x);
            force.y                -= (c21.y + c23.y);
            force.z                -= (c21.z + c23.z);
            cSim.pForce4[offset]    = force;
            offset                  = atom1.z + atom2.y * cSim.stride;
233
            force                   = cSim.pForce4[offset];
Peter Eastman's avatar
Peter Eastman committed
234
235
236
237
238
239
240
            force.x                += c23.x;
            force.y                += c23.y;
            force.z                += c23.z;
            cSim.pForce4[offset]    = force;
        }
        pos += blockDim.x * gridDim.x;
    }
241

Peter Eastman's avatar
Peter Eastman committed
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
    while (pos < cSim.dihedral_offset)
    {
        unsigned int pos1 = pos - cSim.bond_angle_offset;
        if (pos1 < cSim.dihedrals)
        {
            int4   atom1        = cSim.pDihedralID1[pos1];  
            float4 atomA        = cSim.pPosq[atom1.x];
            float4 atomB        = cSim.pPosq[atom1.y];
            float4 atomC        = cSim.pPosq[atom1.z];
            float4 atomD        = cSim.pPosq[atom1.w];            
            A->v0.x             = atomA.x - atomB.x;
            A->v0.y             = atomA.y - atomB.y;
            A->v0.z             = atomA.z - atomB.z;
            A->v1.x             = atomC.x - atomB.x;
            A->v1.y             = atomC.y - atomB.y;
            A->v1.z             = atomC.z - atomB.z;
            A->v2.x             = atomC.x - atomD.x;
            A->v2.y             = atomC.y - atomD.y;
            A->v2.z             = atomC.z - atomD.z; 
            float3 cp0, cp1;
            float dihedralAngle;
            GETDIHEDRALANGLEBETWEENTHREEVECTORS(A->v0, A->v1, A->v2, A->v0, cp0, cp1, dihedralAngle);
            float4 dihedral         = cSim.pDihedralParameter[pos1];
265
            float deltaAngle        = dihedral.z * dihedralAngle - (dihedral.y * LOCAL_HACK_PI / 180.0f);
266
267
268
269
270
271
272
273

	    // ATTENTION: This section leads to a divergent deltaAngle values wrt
	    // forces and energies. We separate the case dihedral.z = n = 0, which
	    // is treated by the calculation of energies via a harmonic potential
/* E */     if (dihedral.z) energy += dihedral.x * (1.0f + cos(deltaAngle));
/* E */     else
	    {
		float deltaAngle    = dihedralAngle - dihedral.y;
274
275
		if (deltaAngle < -LOCAL_HACK_PI) deltaAngle += 2.0f * LOCAL_HACK_PI;
		else if (deltaAngle > LOCAL_HACK_PI) deltaAngle -= 2.0f * LOCAL_HACK_PI;
276
277
278
                energy             += dihedral.x * deltaAngle * deltaAngle;
	    }

Peter Eastman's avatar
Peter Eastman committed
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
            float sinDeltaAngle     = sin(deltaAngle);
            float dEdAngle          = -dihedral.x * dihedral.z * sinDeltaAngle;
            float normCross1        = DOT3(cp0, cp0);
            float normBC            = sqrt(DOT3(A->v1, A->v1));
            float4 ff;
            ff.x                    = (-dEdAngle * normBC) / normCross1;
            float normCross2        = DOT3(cp1, cp1);
            ff.w                    = (dEdAngle * normBC) / normCross2;
            float dp                = 1.0f / DOT3(A->v1, A->v1);
            ff.y                    = DOT3(A->v0, A->v1) * dp;
            ff.z                    = DOT3(A->v2, A->v1) * dp;
            int4  atom2             = cSim.pDihedralID2[pos1];   
            float3 internalF0;
            float3 internalF3;
            float3 s;
            
//            printf("%4d: %9.4f %9.4f %9.4f %9.4f\n", pos1, ff.x, ff.y, ff.z, ff.w);  
            unsigned int offset                 = atom1.x + atom2.x * cSim.stride;
297
            float4 force                        = cSim.pForce4[offset]; 
Peter Eastman's avatar
Peter Eastman committed
298
299
300
301
302
303
304
305
306
307
            internalF0.x                        = ff.x * cp0.x; 
            force.x                            += internalF0.x;
            internalF0.y                        = ff.x * cp0.y;
            force.y                            += internalF0.y;
            internalF0.z                        = ff.x * cp0.z;       
            force.z                            += internalF0.z;
            cSim.pForce4[offset]                = force;
            
            //printf("%4d - 0: %9.4f %9.4f %9.4f\n", pos1, cSim.pForce[offset], cSim.pForce[offset + cSim.stride], cSim.pForce[offset + cSim.stride2]);
            offset                              = atom1.w + atom2.w * cSim.stride;
308
            force                               = cSim.pForce4[offset];
Peter Eastman's avatar
Peter Eastman committed
309
310
311
312
313
314
315
316
317
318
319
320
321
            internalF3.x                        = ff.w * cp1.x;
            force.x                            += internalF3.x;
            internalF3.y                        = ff.w * cp1.y;
            force.y                            += internalF3.y;
            internalF3.z                        = ff.w * cp1.z;
            force.z                            += internalF3.z;
            cSim.pForce4[offset]                = force;
            
           // printf("%4d - 3: %9.4f %9.4f %9.4f\n", pos1, cSim.pForce[offset], cSim.pForce[offset + cSim.stride], cSim.pForce[offset + cSim.stride2]);
            s.x                                 = ff.y * internalF0.x - ff.z * internalF3.x;   
            s.y                                 = ff.y * internalF0.y - ff.z * internalF3.y;  
            s.z                                 = ff.y * internalF0.z - ff.z * internalF3.z;        
            offset                              = atom1.y + atom2.y * cSim.stride;
322
            force                               = cSim.pForce4[offset];
Peter Eastman's avatar
Peter Eastman committed
323
324
325
326
327
328
329
            force.x                            += -internalF0.x + s.x;
            force.y                            += -internalF0.y + s.y;
            force.z                            += -internalF0.z + s.z;
            cSim.pForce4[offset]                = force;
            
            //printf("%4d - 1: %9.4f %9.4f %9.4f\n", pos1, cSim.pForce[offset], cSim.pForce[offset + cSim.stride], cSim.pForce[offset + cSim.stride2]);
            offset                              = atom1.z + atom2.z * cSim.stride;
330
            force                               = cSim.pForce4[offset];
Peter Eastman's avatar
Peter Eastman committed
331
332
333
334
335
            force.x                            += -internalF3.x - s.x;
            force.y                            += -internalF3.y - s.y;
            force.z                            += -internalF3.z - s.z;
            cSim.pForce4[offset]                = force;
            //printf("%4d - 2: %9.4f %9.4f %9.4f\n", pos1, cSim.pForce[offset], cSim.pForce[offset + cSim.stride], cSim.pForce[offset + cSim.stride2]);
336
        }
Peter Eastman's avatar
Peter Eastman committed
337
338
339
        pos += blockDim.x * gridDim.x;
    }

340
    // Ryckaert Bellemans dihedrals
Peter Eastman's avatar
Peter Eastman committed
341
342
343
344
345
    while (pos < cSim.rb_dihedral_offset)
    {
        unsigned int pos1 = pos - cSim.dihedral_offset;
        if (pos1 < cSim.rb_dihedrals)
        {
346
            int4   atom1        = cSim.pRbDihedralID1[pos1];
Peter Eastman's avatar
Peter Eastman committed
347
348
349
            float4 atomA        = cSim.pPosq[atom1.x];
            float4 atomB        = cSim.pPosq[atom1.y];
            float4 atomC        = cSim.pPosq[atom1.z];
350
            float4 atomD        = cSim.pPosq[atom1.w];
Peter Eastman's avatar
Peter Eastman committed
351
352
353
354
355
356
357
358
            A->v0.x             = atomA.x - atomB.x;
            A->v0.y             = atomA.y - atomB.y;
            A->v0.z             = atomA.z - atomB.z;
            A->v1.x             = atomC.x - atomB.x;
            A->v1.y             = atomC.y - atomB.y;
            A->v1.z             = atomC.z - atomB.z;
            A->v2.x             = atomC.x - atomD.x;
            A->v2.y             = atomC.y - atomD.y;
359
            A->v2.z             = atomC.z - atomD.z;
Peter Eastman's avatar
Peter Eastman committed
360
361
            float3 cp0, cp1;
            float dihedralAngle, cosPhi;
362
363
364
      //      printf("%4d - 0 : %9.4f %9.4f %9.4f\n", pos1, A->v0.x, A->v0.y, A->v0.z);
      //      printf("%4d - 1 : %9.4f %9.4f %9.4f\n", pos1, A->v1.x, A->v1.y, A->v1.z);
      //      printf("%4d - 2 : %9.4f %9.4f %9.4f\n", pos1, A->v2.x, A->v2.y, A->v2.z);
Peter Eastman's avatar
Peter Eastman committed
365
366
367
            GETDIHEDRALANGLECOSINEBETWEENTHREEVECTORS(A->v0, A->v1, A->v2, A->v0, cp0, cp1, dihedralAngle, cosPhi);
            if (dihedralAngle < 0.0f )
            {
368
                dihedralAngle += LOCAL_HACK_PI;
369
370
            }
            else
Peter Eastman's avatar
Peter Eastman committed
371
            {
372
                dihedralAngle -= LOCAL_HACK_PI;
Peter Eastman's avatar
Peter Eastman committed
373
374
375
376
377
378
379
            }
            cosPhi                  = -cosPhi;
         //   printf("%4d: %9.4f %9.4f\n", pos1, dihedralAngle, cosPhi);
            float4 dihedral1        = cSim.pRbDihedralParameter1[pos1];
            float2 dihedral2        = cSim.pRbDihedralParameter2[pos1];
            float cosFactor         = cosPhi;
            float dEdAngle          = -dihedral1.y;
380
381
382

/* E */     float rb_energy         = dihedral1.x;
            rb_energy              += dihedral1.y * cosFactor;
Peter Eastman's avatar
Peter Eastman committed
383
384
385
386
387
        //    printf("%4d - 1: %9.4f %9.4f\n", pos1, dEdAngle, 1.0f);
            dEdAngle               -= 2.0f * dihedral1.z * cosFactor;
       //     printf("%4d - 2: %9.4f %9.4f\n", pos1, dEdAngle, cosFactor);
            cosFactor              *= cosPhi;
            dEdAngle               -= 3.0f * dihedral1.w * cosFactor;
388
389
            rb_energy              += dihedral1.z * cosFactor;
    //       printf("%4d - 3: %9.4f %9.4f\n", pos1, dEdAngle, cosFactor);
Peter Eastman's avatar
Peter Eastman committed
390
391
            cosFactor              *= cosPhi;
            dEdAngle               -= 4.0f * dihedral2.x * cosFactor;
392
393
            rb_energy              += dihedral1.w * cosFactor;
  //         printf("%4d - 4: %9.4f %9.4f\n", pos1, dEdAngle, cosFactor);
Peter Eastman's avatar
Peter Eastman committed
394
395
            cosFactor              *= cosPhi;
            dEdAngle               -= 5.0f * dihedral2.y * cosFactor;
396
397
398
            rb_energy              += dihedral2.x * cosFactor;
            rb_energy              += dihedral2.y * cosFactor * cosPhi;
/* E */     energy                 += rb_energy;
Peter Eastman's avatar
Peter Eastman committed
399
 //           printf("%4d - 5: %9.4f %9.4f\n", pos1, dEdAngle, cosFactor);
400
            dEdAngle               *= sin(dihedralAngle);
Peter Eastman's avatar
Peter Eastman committed
401
//            printf("%4d - f: %9.4f\n", pos1, dEdAngle);
402

Peter Eastman's avatar
Peter Eastman committed
403
404
405
406
407
408
409
410
411
            float normCross1        = DOT3(cp0, cp0);
            float normBC            = sqrt(DOT3(A->v1, A->v1));
            float4 ff;
            ff.x                    = (-dEdAngle * normBC) / normCross1;
            float normCross2        = DOT3(cp1, cp1);
            ff.w                    = (dEdAngle * normBC) / normCross2;
            float dp                = 1.0f / DOT3(A->v1, A->v1);
            ff.y                    = DOT3(A->v0, A->v1) * dp;
            ff.z                    = DOT3(A->v2, A->v1) * dp;
412
            int4  atom2             = cSim.pRbDihedralID2[pos1];
Peter Eastman's avatar
Peter Eastman committed
413
414
415
            float3 internalF0;
            float3 internalF3;
            float3 s;
416
417

//            printf("%4d: %9.4f %9.4f %9.4f %9.4f\n", pos1, ff.x, ff.y, ff.z, ff.w);
Peter Eastman's avatar
Peter Eastman committed
418
            unsigned int offset                 = atom1.x + atom2.x * cSim.stride;
419
420
            float4 force                        = cSim.pForce4[offset];
            internalF0.x                        = ff.x * cp0.x;
Peter Eastman's avatar
Peter Eastman committed
421
422
423
            force.x                            += internalF0.x;
            internalF0.y                        = ff.x * cp0.y;
            force.y                            += internalF0.y;
424
            internalF0.z                        = ff.x * cp0.z;
Peter Eastman's avatar
Peter Eastman committed
425
426
            force.z                            += internalF0.z;
            cSim.pForce4[offset]                = force;
427

Peter Eastman's avatar
Peter Eastman committed
428
429
 //           printf("%4d - 0: %9.4f %9.4f %9.4f\n", pos1, cSim.pForce[offset], cSim.pForce[offset + cSim.stride], cSim.pForce[offset + cSim.stride2]);
            offset                              = atom1.w + atom2.w * cSim.stride;
430
            force                               = cSim.pForce4[offset];
Peter Eastman's avatar
Peter Eastman committed
431
432
433
434
435
436
437
            internalF3.x                        = ff.w * cp1.x;
            force.x                            += internalF3.x;
            internalF3.y                        = ff.w * cp1.y;
            force.y                            += internalF3.y;
            internalF3.z                        = ff.w * cp1.z;
            force.z                            += internalF3.z;
            cSim.pForce4[offset]                = force;
438

Peter Eastman's avatar
Peter Eastman committed
439
   //         printf("%4d - 3: %9.4f %9.4f %9.4f\n", pos1, cSim.pForce[offset], cSim.pForce[offset + cSim.stride], cSim.pForce[offset + cSim.stride2]);
440
441
442
            s.x                                 = ff.y * internalF0.x - ff.z * internalF3.x;
            s.y                                 = ff.y * internalF0.y - ff.z * internalF3.y;
            s.z                                 = ff.y * internalF0.z - ff.z * internalF3.z;
Peter Eastman's avatar
Peter Eastman committed
443
            offset                              = atom1.y + atom2.y * cSim.stride;
444
            force                               = cSim.pForce4[offset];
Peter Eastman's avatar
Peter Eastman committed
445
446
447
448
449
450
            force.x                            += -internalF0.x + s.x;
            force.y                            += -internalF0.y + s.y;
            force.z                            += -internalF0.z + s.z;
            cSim.pForce4[offset]                = force;
     //       printf("%4d - 1: %9.4f %9.4f %9.4f\n", pos1, cSim.pForce[offset], cSim.pForce[offset + cSim.stride], cSim.pForce[offset + cSim.stride2]);
            offset                              = atom1.z + atom2.z * cSim.stride;
451
            force                               = cSim.pForce4[offset];
Peter Eastman's avatar
Peter Eastman committed
452
453
454
455
456
            force.x                            += -internalF3.x - s.x;
            force.y                            += -internalF3.y - s.y;
            force.z                            += -internalF3.z - s.z;
            cSim.pForce4[offset]                = force;
     //       printf("%4d - 2: %9.4f %9.4f %9.4f\n", pos1, cSim.pForce[offset], cSim.pForce[offset + cSim.stride], cSim.pForce[offset + cSim.stride2]);
457
        }         
Peter Eastman's avatar
Peter Eastman committed
458
459
460
        pos += blockDim.x * gridDim.x;
    }   

461
    if (cSim.nonbondedMethod == NO_CUTOFF || cSim.nonbondedMethod == EWALD || cSim.nonbondedMethod == PARTICLE_MESH_EWALD)
462
463
    {
        while (pos < cSim.LJ14_offset)
Peter Eastman's avatar
Peter Eastman committed
464
        {
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
            unsigned int pos1       = pos - cSim.rb_dihedral_offset;
            if (pos1 < cSim.LJ14s)
            {
                int4 atom               = cSim.pLJ14ID[pos1];
                float4 LJ14             = cSim.pLJ14Parameter[pos1];
                float4 a1               = cSim.pPosq[atom.x];
                float4 a2               = cSim.pPosq[atom.y];
                float3 d;
                d.x                     = a1.x - a2.x;
                d.y                     = a1.y - a2.y;
                d.z                     = a1.z - a2.z;
                float r2                = DOT3(d, d);
                float inverseR          = 1.0f / sqrt(r2);
                float sig2              = inverseR * LJ14.y;
                sig2                   *= sig2;
                float sig6              = sig2 * sig2 * sig2;
                float dEdR              = LJ14.x * (12.0f * sig6 - 6.0f) * sig6;
482
                /* E */
483
                energy                 += LJ14.x * (sig6 - 1.0f) * sig6;
484
485
                energy                 += LJ14.z * inverseR;

486
487
488
489
                dEdR                   += LJ14.z * inverseR;
                dEdR                   *= inverseR * inverseR;
                unsigned int offsetA    = atom.x + atom.z * cSim.stride;
                unsigned int offsetB    = atom.y + atom.w * cSim.stride;
490
491
                float4 forceA           = cSim.pForce4[offsetA];
                float4 forceB           = cSim.pForce4[offsetB];
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
                d.x                    *= dEdR;
                d.y                    *= dEdR;
                d.z                    *= dEdR;
                forceA.x               += d.x;
                forceA.y               += d.y;
                forceA.z               += d.z;
                forceB.x               -= d.x;
                forceB.y               -= d.y;
                forceB.z               -= d.z;
                cSim.pForce4[offsetA]   = forceA;
                cSim.pForce4[offsetB]   = forceB;
            }
            pos                    += blockDim.x * gridDim.x;
        }
    }
    else if (cSim.nonbondedMethod == CUTOFF)
    {
509
        float LJ14_energy;
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
        while (pos < cSim.LJ14_offset)
        {
            unsigned int pos1       = pos - cSim.rb_dihedral_offset;
            if (pos1 < cSim.LJ14s)
            {
                int4 atom               = cSim.pLJ14ID[pos1];
                float4 LJ14             = cSim.pLJ14Parameter[pos1];
                float4 a1               = cSim.pPosq[atom.x];
                float4 a2               = cSim.pPosq[atom.y];
                float3 d;
                d.x                     = a1.x - a2.x;
                d.y                     = a1.y - a2.y;
                d.z                     = a1.z - a2.z;
                float r2                = DOT3(d, d);
                float inverseR          = 1.0f / sqrt(r2);
                float sig2              = inverseR * LJ14.y;
                sig2                   *= sig2;
                float sig6              = sig2 * sig2 * sig2;
528
529
530
531
                float dEdR              = LJ14.x * (12.0f * sig6 - 6.0f) * sig6;                
                /* E */
                LJ14_energy             = LJ14.x * (sig6 - 1.0f) * sig6;
                LJ14_energy            += LJ14.z * (inverseR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);
532
533
534
                dEdR                   += LJ14.z * (inverseR - 2.0f * cSim.reactionFieldK * r2);
                dEdR                   *= inverseR * inverseR;
                if (r2 > cSim.nonbondedCutoffSqr)
535
                {                   
536
                    dEdR = 0.0f;
537
538
                    /* E */
                    LJ14_energy = 0.0f;
539
                }
540
541
542
                /* E */
                energy                 += LJ14_energy;
 
543
544
                unsigned int offsetA    = atom.x + atom.z * cSim.stride;
                unsigned int offsetB    = atom.y + atom.w * cSim.stride;
545
546
                float4 forceA           = cSim.pForce4[offsetA];
                float4 forceB           = cSim.pForce4[offsetB];
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
                d.x                    *= dEdR;
                d.y                    *= dEdR;
                d.z                    *= dEdR;
                forceA.x               += d.x;
                forceA.y               += d.y;
                forceA.z               += d.z;
                forceB.x               -= d.x;
                forceB.y               -= d.y;
                forceB.z               -= d.z;
                cSim.pForce4[offsetA]   = forceA;
                cSim.pForce4[offsetB]   = forceB;
            }
            pos                    += blockDim.x * gridDim.x;
        }
    }
    else if (cSim.nonbondedMethod == PERIODIC)
    {
564
        float LJ14_energy;
565
566
567
568
569
570
571
572
573
574
575
576
577
578
        while (pos < cSim.LJ14_offset)
        {
            unsigned int pos1       = pos - cSim.rb_dihedral_offset;
            if (pos1 < cSim.LJ14s)
            {
                int4 atom               = cSim.pLJ14ID[pos1];
                float4 LJ14             = cSim.pLJ14Parameter[pos1];
                float4 a1               = cSim.pPosq[atom.x];
                float4 a2               = cSim.pPosq[atom.y];
                float3 d;
                d.x                     = a1.x - a2.x;
                d.y                     = a1.y - a2.y;
                d.z                     = a1.z - a2.z;
                d.x                     -= floor(d.x/cSim.periodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
579
580
                d.y                     -= floor(d.y/cSim.periodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
                d.z                     -= floor(d.z/cSim.periodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
581
582
583
584
585
586
                float r2                = DOT3(d, d);
                float inverseR          = 1.0f / sqrt(r2);
                float sig2              = inverseR * LJ14.y;
                sig2                   *= sig2;
                float sig6              = sig2 * sig2 * sig2;
                float dEdR              = LJ14.x * (12.0f * sig6 - 6.0f) * sig6;
587
588
589
590
                /* E */
                LJ14_energy             = LJ14.x * (sig6 - 1.0f) * sig6;
                LJ14_energy            += LJ14.z * (inverseR + cSim.reactionFieldK * r2 - cSim.reactionFieldC);

591
592
593
594
595
                dEdR                   += LJ14.z * (inverseR - 2.0f * cSim.reactionFieldK * r2);
                dEdR                   *= inverseR * inverseR;
                if (r2 > cSim.nonbondedCutoffSqr)
                {
                    dEdR = 0.0f;
596
597
                    /* E */
                    LJ14_energy = 0.0f;
598
                }
599
600
601
                /* E */
                energy                 += LJ14_energy;

602
603
                unsigned int offsetA    = atom.x + atom.z * cSim.stride;
                unsigned int offsetB    = atom.y + atom.w * cSim.stride;
604
605
                float4 forceA           = cSim.pForce4[offsetA];
                float4 forceB           = cSim.pForce4[offsetB];
606
607
608
609
610
611
612
613
614
615
616
617
618
619
                d.x                    *= dEdR;
                d.y                    *= dEdR;
                d.z                    *= dEdR;
                forceA.x               += d.x;
                forceA.y               += d.y;
                forceA.z               += d.z;
                forceB.x               -= d.x;
                forceB.y               -= d.y;
                forceB.z               -= d.z;
                cSim.pForce4[offsetA]   = forceA;
                cSim.pForce4[offsetB]   = forceB;
            }
            pos                    += blockDim.x * gridDim.x;
        }
Peter Eastman's avatar
Peter Eastman committed
620
    }
621
    cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += energy;
Peter Eastman's avatar
Peter Eastman committed
622
623
}

624

625
void kCalculateLocalForces(gpuContext gpu)
Peter Eastman's avatar
Peter Eastman committed
626
627
{
  //  printf("kCalculateLocalForces\n");
628
    kCalculateLocalForces_kernel<<<gpu->sim.blocks, gpu->sim.localForces_threads_per_block, gpu->sim.localForces_threads_per_block * sizeof(Vectors)>>>();
Peter Eastman's avatar
Peter Eastman committed
629
630
631
    LAUNCHERROR("kCalculateLocalForces");
}