kCalculateCustomNonbondedForces.cu 14.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
/* -------------------------------------------------------------------------- *
 *                                   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:                                                              *
 *                                                                            *
 * 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.                                        *
 *                                                                            *
 * 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.                        *
 *                                                                            *
 * 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/>.      *
 * -------------------------------------------------------------------------- */

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

#include "gputypes.h"
#include "cudatypes.h"

#define UNROLLXX 0
#define UNROLLXY 0

struct Atom {
    float x;
    float y;
    float z;
    float4 params;
    float fx;
    float fy;
    float fz;
};

static __constant__ cudaGmxSimulation cSim;
static __constant__ Expression<128> forceExp;
static __constant__ Expression<128> energyExp;
static __constant__ Expression<64> combiningRules[4];
56
static __constant__ float globalParams[8];
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92

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

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

void SetCustomNonbondedForceExpression(const Expression<128>& expression)
{
    cudaError_t status;
    status = cudaMemcpyToSymbol(forceExp, &expression, sizeof(forceExp));
    RTERROR(status, "SetCustomNonbondedForceExpression: cudaMemcpyToSymbol failed");
}

void SetCustomNonbondedEnergyExpression(const Expression<128>& expression)
{
    cudaError_t status;
    status = cudaMemcpyToSymbol(energyExp, &expression, sizeof(energyExp));
    RTERROR(status, "SetCustomNonbondedEnergyExpression: cudaMemcpyToSymbol failed");
}

void SetCustomNonbondedCombiningRules(const Expression<64>* expressions)
{
    cudaError_t status;
    status = cudaMemcpyToSymbol(combiningRules, expressions, sizeof(combiningRules));
    RTERROR(status, "SetCustomNonbondedCombiningRules: cudaMemcpyToSymbol failed");
}

93
94
95
96
97
98
99
void SetCustomNonbondedGlobalParams(const vector<float>& paramValues)
{
    cudaError_t status;
    status = cudaMemcpyToSymbol(globalParams, &paramValues[0], sizeof(globalParams));
    RTERROR(status, "SetCustomNonbondedGlobalParams: cudaMemcpyToSymbol failed");
}

100
101
#define STACK(y) stack[(y)*blockDim.x+threadIdx.x]

102
103
104
105
106
107
108
109
110
template<int SIZE>
__device__ float kEvaluateExpression_kernel(Expression<SIZE>* expression, float* stack, float var0, float4 vars1, float4 vars2)
{
    int stackPointer = -1;
    for (int i = 0; i < expression->length; i++)
    {
        switch (expression->op[i])
        {
            case CONSTANT:
111
                STACK(++stackPointer) = expression->arg[i];
112
113
                break;
            case VARIABLE0:
114
                STACK(++stackPointer) = var0;
115
116
                break;
            case VARIABLE1:
117
                STACK(++stackPointer) = vars1.x;
118
                break;
119
            case VARIABLE2:
120
                STACK(++stackPointer) = vars1.y;
121
122
                break;
            case VARIABLE3:
123
                STACK(++stackPointer) = vars1.z;
124
125
                break;
            case VARIABLE4:
126
                STACK(++stackPointer) = vars1.w;
127
                break;
128
            case VARIABLE5:
129
                STACK(++stackPointer) = vars2.x;
130
                break;
131
            case VARIABLE6:
132
                STACK(++stackPointer) = vars2.y;
133
134
                break;
            case VARIABLE7:
135
                STACK(++stackPointer) = vars2.z;
136
137
                break;
            case VARIABLE8:
138
                STACK(++stackPointer) = vars2.w;
139
                break;
140
            case GLOBAL:
141
                STACK(++stackPointer) = globalParams[(int) expression->arg[i]];
142
                break;
143
144
            case ADD:
            {
145
146
                float temp = STACK(stackPointer);
                STACK(--stackPointer) += temp;
147
148
149
150
                break;
            }
            case SUBTRACT:
            {
151
152
                float temp = STACK(stackPointer);
                STACK(stackPointer) = temp-STACK(--stackPointer);
153
154
155
156
                break;
            }
            case MULTIPLY:
            {
157
158
                float temp = STACK(stackPointer);
                STACK(--stackPointer) *= temp;
159
160
161
162
                break;
            }
            case DIVIDE:
            {
163
164
                float temp = STACK(stackPointer);
                STACK(stackPointer) = temp/STACK(--stackPointer);
165
166
167
168
                break;
            }
            case POWER:
            {
169
170
                float temp = STACK(stackPointer);
                STACK(stackPointer) = pow(temp, STACK(--stackPointer));
171
172
173
                break;
            }
            case NEGATE:
174
                STACK(stackPointer) *= -1.0f;
175
176
                break;
            case SQRT:
177
                STACK(stackPointer) = sqrt(STACK(stackPointer));
178
179
                break;
            case EXP:
180
                STACK(stackPointer) = exp(STACK(stackPointer));
181
182
                break;
            case LOG:
183
                STACK(stackPointer) = log(STACK(stackPointer));
184
185
                break;
            case SIN:
186
                STACK(stackPointer) = sin(STACK(stackPointer));
187
188
                break;
            case COS:
189
                STACK(stackPointer) = cos(STACK(stackPointer));
190
191
                break;
            case SEC:
192
                STACK(stackPointer) = 1.0f/cos(STACK(stackPointer));
193
194
                break;
            case CSC:
195
                STACK(stackPointer) = 1.0f/sin(STACK(stackPointer));
196
197
                break;
            case TAN:
198
                STACK(stackPointer) = tan(STACK(stackPointer));
199
200
                break;
            case COT:
201
                STACK(stackPointer) = 1.0f/tan(STACK(stackPointer));
202
203
                break;
            case ASIN:
204
                STACK(stackPointer) = asin(STACK(stackPointer));
205
206
                break;
            case ACOS:
207
                STACK(stackPointer) = acos(STACK(stackPointer));
208
209
                break;
            case ATAN:
210
                STACK(stackPointer) = atan(STACK(stackPointer));
211
212
213
                break;
            case SQUARE:
            {
214
215
                float temp = STACK(stackPointer);
                STACK(stackPointer) *= temp;
216
217
218
219
                break;
            }
            case CUBE:
            {
220
221
                float temp = STACK(stackPointer);
                STACK(stackPointer) *= temp*temp;
222
223
224
                break;
            }
            case RECIPROCAL:
225
                STACK(stackPointer) = 1.0f/STACK(stackPointer);
226
227
                break;
            case INCREMENT:
228
                STACK(stackPointer) += 1.0f;
229
230
                break;
            case DECREMENT:
231
                STACK(stackPointer) -= 1.0f;
232
233
234
                break;
        }
    }
235
    return STACK(stackPointer);
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
}

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

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

// Include versions of the kernels with cutoffs.

#undef METHOD_NAME
#undef USE_OUTPUT_BUFFER_PER_WARP
#define USE_CUTOFF
#define METHOD_NAME(a, b) a##Cutoff##b
#include "kCalculateCustomNonbondedForces.h"
#define USE_OUTPUT_BUFFER_PER_WARP
#undef METHOD_NAME
#define METHOD_NAME(a, b) a##CutoffByWarp##b
#include "kCalculateCustomNonbondedForces.h"

// Include versions of the kernels with periodic boundary conditions.

#undef METHOD_NAME
#undef USE_OUTPUT_BUFFER_PER_WARP
#define USE_PERIODIC
#define METHOD_NAME(a, b) a##Periodic##b
#include "kCalculateCustomNonbondedForces.h"
#define USE_OUTPUT_BUFFER_PER_WARP
#undef METHOD_NAME
#define METHOD_NAME(a, b) a##PeriodicByWarp##b
#include "kCalculateCustomNonbondedForces.h"

__global__ void kFindBlockBoundsCutoff_kernel();
__global__ void kFindBlocksWithInteractionsCutoff_kernel();
__global__ void kFindInteractionsWithinBlocksCutoff_kernel(unsigned int* workUnit);
__global__ void kFindBlockBoundsPeriodic_kernel();
__global__ void kFindBlocksWithInteractionsPeriodic_kernel();
__global__ void kFindInteractionsWithinBlocksPeriodic_kernel(unsigned int* workUnit);

void kCalculateCustomNonbondedForces(gpuContext gpu, bool neighborListValid)
{
//    printf("kCalculateCustomNonbondedCutoffForces\n");
    CUDPPResult result;
282
283
284
285
286
287
288
    int sharedPerThread = sizeof(Atom)+gpu->sim.customExpressionStackSize*sizeof(float);
    if (gpu->sim.customNonbondedMethod != NO_CUTOFF)
        sharedPerThread += sizeof(float3);
    int threads = gpu->sim.nonbond_threads_per_block;
    int maxThreads = 16380/sharedPerThread;
    if (threads > maxThreads)
        threads = (maxThreads/32)*32;
289
290
291
292
    switch (gpu->sim.customNonbondedMethod)
    {
        case NO_CUTOFF:
            if (gpu->bOutputBufferPerWarp)
293
                kCalculateCustomNonbondedN2ByWarpForces_kernel<<<gpu->sim.nonbond_blocks, threads, sharedPerThread*threads>>>(gpu->sim.pWorkUnit);
294
            else
295
                kCalculateCustomNonbondedN2Forces_kernel<<<gpu->sim.nonbond_blocks, threads, sharedPerThread*threads>>>(gpu->sim.pWorkUnit);
296
            LAUNCHERROR("kCalculateCustomNonbondedN2Forces");
297
298
299
            kCalculateCustomNonbondedN2Exceptions_kernel<<<gpu->sim.blocks, gpu->sim.custom_exception_threads_per_block,
                    gpu->sim.customExpressionStackSize*sizeof(float)*gpu->sim.custom_exception_threads_per_block>>>();
            LAUNCHERROR("kCalculateCustomNonbondedN2Exceptions");
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
            break;
        case CUTOFF:
            if (!neighborListValid)
            {
                kFindBlockBoundsCutoff_kernel<<<(gpu->psGridBoundingBox->_length+63)/64, 64>>>();
                LAUNCHERROR("kFindBlockBoundsCutoff");
                kFindBlocksWithInteractionsCutoff_kernel<<<gpu->sim.interaction_blocks, gpu->sim.interaction_threads_per_block>>>();
                LAUNCHERROR("kFindBlocksWithInteractionsCutoff");
                result = cudppCompact(gpu->cudpp, gpu->sim.pInteractingWorkUnit, gpu->sim.pInteractionCount,
                        gpu->sim.pWorkUnit, gpu->sim.pInteractionFlag, gpu->sim.workUnits);
                if (result != CUDPP_SUCCESS)
                {
                    printf("Error in cudppCompact: %d\n", result);
                    exit(-1);
                }
                kFindInteractionsWithinBlocksCutoff_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
                        sizeof(unsigned int)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
            }
            if (gpu->bOutputBufferPerWarp)
319
                kCalculateCustomNonbondedCutoffByWarpForces_kernel<<<gpu->sim.nonbond_blocks, threads, sharedPerThread*threads>>>(gpu->sim.pInteractingWorkUnit);
320
            else
321
                kCalculateCustomNonbondedCutoffForces_kernel<<<gpu->sim.nonbond_blocks, threads, sharedPerThread*threads>>>(gpu->sim.pInteractingWorkUnit);
322
            LAUNCHERROR("kCalculateCustomNonbondedCutoffForces");
323
324
325
            kCalculateCustomNonbondedCutoffExceptions_kernel<<<gpu->sim.blocks, gpu->sim.custom_exception_threads_per_block,
                    gpu->sim.customExpressionStackSize*sizeof(float)*gpu->sim.custom_exception_threads_per_block>>>();
            LAUNCHERROR("kCalculateCustomNonbondedCutoffExceptions");
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
            break;
        case PERIODIC:
            if (!neighborListValid)
            {
                kFindBlockBoundsPeriodic_kernel<<<(gpu->psGridBoundingBox->_length+63)/64, 64>>>();
                LAUNCHERROR("kFindBlockBoundsPeriodic");
                kFindBlocksWithInteractionsPeriodic_kernel<<<gpu->sim.interaction_blocks, gpu->sim.interaction_threads_per_block>>>();
                LAUNCHERROR("kFindBlocksWithInteractionsPeriodic");
                result = cudppCompact(gpu->cudpp, gpu->sim.pInteractingWorkUnit, gpu->sim.pInteractionCount,
                        gpu->sim.pWorkUnit, gpu->sim.pInteractionFlag, gpu->sim.workUnits);
                if (result != CUDPP_SUCCESS)
                {
                    printf("Error in cudppCompact: %d\n", result);
                    exit(-1);
                }
                kFindInteractionsWithinBlocksPeriodic_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
                        sizeof(unsigned int)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
            }
            if (gpu->bOutputBufferPerWarp)
345
                kCalculateCustomNonbondedPeriodicByWarpForces_kernel<<<gpu->sim.nonbond_blocks, threads, sharedPerThread*threads>>>(gpu->sim.pInteractingWorkUnit);
346
            else
347
                kCalculateCustomNonbondedPeriodicForces_kernel<<<gpu->sim.nonbond_blocks, threads, sharedPerThread*threads>>>(gpu->sim.pInteractingWorkUnit);
348
            LAUNCHERROR("kCalculateCustomNonbondedPeriodicForces");
349
350
351
            kCalculateCustomNonbondedPeriodicExceptions_kernel<<<gpu->sim.blocks, gpu->sim.custom_exception_threads_per_block,
                    gpu->sim.customExpressionStackSize*sizeof(float)*gpu->sim.custom_exception_threads_per_block>>>();
            LAUNCHERROR("kCalculateCustomNonbondedPeriodicExceptions");
352
353
354
            break;
    }
}