CudaKernels.cpp 32.2 KB
Newer Older
1
2
3
4
5
6
7
8
/* -------------------------------------------------------------------------- *
 *                                   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.               *
 *                                                                            *
9
 * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
10
11
12
 * Authors: 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.                                        *
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.                        *
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/>.      *
25
26
27
28
 * -------------------------------------------------------------------------- */

#include "CudaKernels.h"
#include "CudaStreamImpl.h"
29
#include "openmm/LangevinIntegrator.h"
30
#include "openmm/Context.h"
31
#include "ReferencePlatform.h"
32
#include "openmm/internal/ContextImpl.h"
33
#include "kernels/gputypes.h"
34
#include "kernels/cudaKernels.h"
35
36
37
38
39
40
41
#include <cmath>

extern "C" int gpuSetConstants( gpuContext gpu );

using namespace OpenMM;
using namespace std;

42
static void calcForces(ContextImpl& context, CudaPlatform::PlatformData& data) {
43
    _gpuContext* gpu = data.gpu;
44
    if (data.nonbondedMethod != NO_CUTOFF && data.computeForceCount%100 == 0)
45
        gpuReorderAtoms(gpu);
46
    data.computeForceCount++;
47
48
    kClearForces(gpu);
    if (gpu->bIncludeGBSA) {
49
        gpu->bRecalculateBornRadii = true;
50
51
52
53
        kCalculateCDLJObcGbsaForces1(gpu);
        kReduceObcGbsaBornForces(gpu);
        kCalculateObcGbsaForces2(gpu);
    }
54
    else if (data.hasNonbonded)
55
        kCalculateCDLJForces(gpu);
56
57
    if (data.hasCustomNonbonded)
        kCalculateCustomNonbondedForces(gpu, data.hasNonbonded);
58
    kCalculateLocalForces(gpu);
59
    kReduceForces(gpu);
60
61
}

62
63
static double calcEnergy(ContextImpl& context, CudaPlatform::PlatformData& data, System& system) {
    _gpuContext* gpu = data.gpu;
64
65
66
67
68
69
70
71
72
    if (data.nonbondedMethod != NO_CUTOFF && data.stepCount%100 == 0)
        gpuReorderAtoms(gpu);
    data.stepCount++;
    kClearEnergy(gpu);
    if (gpu->bIncludeGBSA) {
        gpu->bRecalculateBornRadii = true;
        kCalculateCDLJObcGbsaForces1(gpu);
        kReduceObcGbsaBornForces(gpu);
        kCalculateObcGbsaForces2(gpu);
73
    }
74
75
76
77
78
79
80
81
    else if (data.hasNonbonded)
        kCalculateCDLJForces(gpu);
    if (data.hasCustomNonbonded)
        kCalculateCustomNonbondedForces(gpu, data.hasNonbonded);
    kCalculateLocalForces(gpu);
    if (gpu->bIncludeGBSA)
        kReduceBornSumAndForces(gpu);
    return kReduceEnergy(gpu)+data.ewaldSelfEnergy;
82
83
}

84
85
86
void CudaInitializeForcesKernel::initialize(const System& system) {
}

87
void CudaInitializeForcesKernel::execute(ContextImpl& context) {
88
89
}

90
91
92
void CudaUpdateTimeKernel::initialize(const System& system) {
}

93
double CudaUpdateTimeKernel::getTime(const ContextImpl& context) const {
94
95
96
    return data.time;
}

97
void CudaUpdateTimeKernel::setTime(ContextImpl& context, double time) {
98
99
100
    data.time = time;
}

101
102
103
104
105
106
107
CudaCalcHarmonicBondForceKernel::~CudaCalcHarmonicBondForceKernel() {
}

void CudaCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasBonds = true;
108
    numBonds = force.getNumBonds();
Peter Eastman's avatar
Peter Eastman committed
109
110
    vector<int> particle1(numBonds);
    vector<int> particle2(numBonds);
111
112
113
114
    vector<float> length(numBonds);
    vector<float> k(numBonds);
    for (int i = 0; i < numBonds; i++) {
        double lengthValue, kValue;
Peter Eastman's avatar
Peter Eastman committed
115
        force.getBondParameters(i, particle1[i], particle2[i], lengthValue, kValue);
116
117
118
        length[i] = (float) lengthValue;
        k[i] = (float) kValue;
    }
Peter Eastman's avatar
Peter Eastman committed
119
    gpuSetBondParameters(data.gpu, particle1, particle2, length, k);
120
121
}

122
void CudaCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
123
124
125
126
    if (data.primaryKernel == this)
        calcForces(context, data);
}

127
double CudaCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
128
    if (data.primaryKernel == this)
129
        return calcEnergy(context, data, system);
130
131
132
133
134
135
136
137
138
139
    return 0.0;
}

CudaCalcHarmonicAngleForceKernel::~CudaCalcHarmonicAngleForceKernel() {
}

void CudaCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasAngles = true;
140
    numAngles = force.getNumAngles();
141
    const float RadiansToDegrees = (float) (180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
142
143
144
    vector<int> particle1(numAngles);
    vector<int> particle2(numAngles);
    vector<int> particle3(numAngles);
145
146
147
148
    vector<float> angle(numAngles);
    vector<float> k(numAngles);
    for (int i = 0; i < numAngles; i++) {
        double angleValue, kValue;
Peter Eastman's avatar
Peter Eastman committed
149
        force.getAngleParameters(i, particle1[i], particle2[i], particle3[i], angleValue, kValue);
150
151
        angle[i] = (float) (angleValue*RadiansToDegrees);
        k[i] = (float) kValue;
152
    }
Peter Eastman's avatar
Peter Eastman committed
153
    gpuSetBondAngleParameters(data.gpu, particle1, particle2, particle3, angle, k);
154
}
155

156
void CudaCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
157
158
159
160
    if (data.primaryKernel == this)
        calcForces(context, data);
}

161
double CudaCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
162
    if (data.primaryKernel == this)
163
        return calcEnergy(context, data, system);
164
165
166
167
168
169
170
171
172
173
174
    return 0.0;
}

CudaCalcPeriodicTorsionForceKernel::~CudaCalcPeriodicTorsionForceKernel() {
}

void CudaCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasPeriodicTorsions = true;
    numTorsions = force.getNumTorsions();
175
    const float RadiansToDegrees = (float)(180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
176
177
178
179
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
180
181
182
183
184
    vector<float> k(numTorsions);
    vector<float> phase(numTorsions);
    vector<int> periodicity(numTorsions);
    for (int i = 0; i < numTorsions; i++) {
        double kValue, phaseValue;
Peter Eastman's avatar
Peter Eastman committed
185
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], periodicity[i], phaseValue, kValue);
186
187
        k[i] = (float) kValue;
        phase[i] = (float) (phaseValue*RadiansToDegrees);
188
    }
Peter Eastman's avatar
Peter Eastman committed
189
    gpuSetDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, k, phase, periodicity);
190
191
}

192
void CudaCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
193
194
195
196
    if (data.primaryKernel == this)
        calcForces(context, data);
}

197
double CudaCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
198
    if (data.primaryKernel == this)
199
        return calcEnergy(context, data, system);
200
201
202
203
204
205
206
207
208
209
210
    return 0.0;
}

CudaCalcRBTorsionForceKernel::~CudaCalcRBTorsionForceKernel() {
}

void CudaCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasRB = true;
    numTorsions = force.getNumTorsions();
Peter Eastman's avatar
Peter Eastman committed
211
212
213
214
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
215
216
217
218
219
220
221
222
    vector<float> c0(numTorsions);
    vector<float> c1(numTorsions);
    vector<float> c2(numTorsions);
    vector<float> c3(numTorsions);
    vector<float> c4(numTorsions);
    vector<float> c5(numTorsions);
    for (int i = 0; i < numTorsions; i++) {
        double c[6];
Peter Eastman's avatar
Peter Eastman committed
223
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], c[0], c[1], c[2], c[3], c[4], c[5]);
224
225
226
227
228
229
        c0[i] = (float) c[0];
        c1[i] = (float) c[1];
        c2[i] = (float) c[2];
        c3[i] = (float) c[3];
        c4[i] = (float) c[4];
        c5[i] = (float) c[5];
230
    }
Peter Eastman's avatar
Peter Eastman committed
231
    gpuSetRbDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
232
233
}

234
void CudaCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
235
236
237
238
    if (data.primaryKernel == this)
        calcForces(context, data);
}

239
double CudaCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
240
    if (data.primaryKernel == this)
241
        return calcEnergy(context, data, system);
242
243
244
245
246
247
    return 0.0;
}

CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() {
}

248
void CudaCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
249
250
251
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasNonbonded = true;
Peter Eastman's avatar
Peter Eastman committed
252
    numParticles = force.getNumParticles();
253
    _gpuContext* gpu = data.gpu;
254
255
256
257

    // Identify which exceptions are 1-4 interactions.

    vector<pair<int, int> > exclusions;
258
    vector<int> exceptions;
259
260
261
262
263
264
    for (int i = 0; i < force.getNumExceptions(); i++) {
        int particle1, particle2;
        double chargeProd, sigma, epsilon;
        force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
        exclusions.push_back(pair<int, int>(particle1, particle2));
        if (chargeProd != 0.0 || epsilon != 0.0)
265
            exceptions.push_back(i);
266
267
    }

268
269
    // Initialize nonbonded interactions.
    
270
    {
Peter Eastman's avatar
Peter Eastman committed
271
272
273
274
        vector<int> particle(numParticles);
        vector<float> c6(numParticles);
        vector<float> c12(numParticles);
        vector<float> q(numParticles);
275
        vector<char> symbol;
Peter Eastman's avatar
Peter Eastman committed
276
277
        vector<vector<int> > exclusionList(numParticles);
        for (int i = 0; i < numParticles; i++) {
278
            double charge, radius, depth;
Peter Eastman's avatar
Peter Eastman committed
279
280
            force.getParticleParameters(i, charge, radius, depth);
            particle[i] = i;
281
282
283
            q[i] = (float) charge;
            c6[i] = (float) (4*depth*pow(radius, 6.0));
            c12[i] = (float) (4*depth*pow(radius, 12.0));
284
285
            exclusionList[i].push_back(i);
        }
286
        for (int i = 0; i < (int)exclusions.size(); i++) {
287
288
289
            exclusionList[exclusions[i].first].push_back(exclusions[i].second);
            exclusionList[exclusions[i].second].push_back(exclusions[i].first);
        }
290
291
        CudaNonbondedMethod method = NO_CUTOFF;
        if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
292
            gpuSetNonbondedCutoff(gpu, (float)force.getCutoffDistance(), force.getReactionFieldDielectric());
293
294
295
296
297
            method = CUTOFF;
        }
        if (force.getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
            Vec3 boxVectors[3];
            force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
298
            gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
299
300
            method = PERIODIC;
        }
301
        if (force.getNonbondedMethod() == NonbondedForce::Ewald || force.getNonbondedMethod() == NonbondedForce::PME) {
302
303
            Vec3 boxVectors[3];
            force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
304
            gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
305
306
307
308
309
310
            double ewaldErrorTol = force.getEwaldErrorTolerance();
            double alpha = (1.0/force.getCutoffDistance())*std::sqrt(-std::log(ewaldErrorTol));
            double mx = boxVectors[0][0]/force.getCutoffDistance();
            double my = boxVectors[1][1]/force.getCutoffDistance();
            double mz = boxVectors[2][2]/force.getCutoffDistance();
            double pi = 3.1415926535897932385;
Peter Eastman's avatar
Peter Eastman committed
311
312
313
            int kmaxx = (int)std::ceil(-(mx/pi)*std::log(ewaldErrorTol));
            int kmaxy = (int)std::ceil(-(my/pi)*std::log(ewaldErrorTol));
            int kmaxz = (int)std::ceil(-(mz/pi)*std::log(ewaldErrorTol));
314
315
316
317
318
319
320
321
322
323
324
            if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
                if (kmaxx%2 == 0)
                    kmaxx++;
                if (kmaxy%2 == 0)
                    kmaxy++;
                if (kmaxz%2 == 0)
                    kmaxz++;
                gpuSetEwaldParameters(gpu, (float) alpha, kmaxx, kmaxy, kmaxz);
                method = EWALD;
            }
            else {
325
326
327
                int gridSizeX = -0.5*kmaxx*std::log(ewaldErrorTol);
                int gridSizeY = -0.5*kmaxy*std::log(ewaldErrorTol);
                int gridSizeZ = -0.5*kmaxz*std::log(ewaldErrorTol);
Peter Eastman's avatar
Peter Eastman committed
328
                gpuSetPMEParameters(gpu, (float) alpha, gridSizeX, gridSizeY, gridSizeZ);
329
330
                method = PARTICLE_MESH_EWALD;
            }
331
        }
332
333
        data.nonbondedMethod = method;
        gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method);
334
335
336
337

        // Compute the Ewald self energy.

        data.ewaldSelfEnergy = 0.0;
338
339
340
341
342
        if (force.getNonbondedMethod() == NonbondedForce::Ewald || force.getNonbondedMethod() == NonbondedForce::PME) {
            double selfEnergyScale = gpu->sim.epsfac*gpu->sim.alphaEwald/std::sqrt(PI);
                for (int i = 0; i < numParticles; i++)
                    data.ewaldSelfEnergy -= selfEnergyScale*q[i]*q[i];
        }
343
344
345
346
    }

    // Initialize 1-4 nonbonded interactions.
    
347
    {
348
349
350
351
352
353
354
355
        int numExceptions = exceptions.size();
        vector<int> particle1(numExceptions);
        vector<int> particle2(numExceptions);
        vector<float> c6(numExceptions);
        vector<float> c12(numExceptions);
        vector<float> q1(numExceptions);
        vector<float> q2(numExceptions);
        for (int i = 0; i < numExceptions; i++) {
356
            double charge, sig, eps;
357
            force.getExceptionParameters(exceptions[i], particle1[i], particle2[i], charge, sig, eps);
358
359
            c6[i] = (float) (4*eps*pow(sig, 6.0));
            c12[i] = (float) (4*eps*pow(sig, 12.0));
Peter Eastman's avatar
Peter Eastman committed
360
361
            q1[i] = (float) charge;
            q2[i] = 1.0f;
362
        }
Peter Eastman's avatar
Peter Eastman committed
363
        gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, particle1, particle2, c6, c12, q1, q2);
364
365
366
    }
}

367
void CudaCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
368
369
    if (data.primaryKernel == this)
        calcForces(context, data);
370
371
}

372
double CudaCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
373
    if (data.primaryKernel == this)
374
        return calcEnergy(context, data, system);
375
    return 0.0;
376
377
}

378
379
380
381
CudaCalcCustomNonbondedForceKernel::~CudaCalcCustomNonbondedForceKernel() {
}

void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
382
    data.primaryKernel = this; // This must always be the primary kernel so it can update the global parameters
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
    data.hasCustomNonbonded = true;
    numParticles = force.getNumParticles();
    _gpuContext* gpu = data.gpu;

    // Identify which exceptions are actual interactions.

    vector<pair<int, int> > exclusions;
    vector<int> exceptions;
    {
        vector<double> parameters;
        for (int i = 0; i < force.getNumExceptions(); i++) {
            int particle1, particle2;
            force.getExceptionParameters(i, particle1, particle2, parameters);
            exclusions.push_back(pair<int, int>(particle1, particle2));
            if (parameters.size() > 0)
                exceptions.push_back(i);
        }
    }

    // Initialize nonbonded interactions.

    vector<int> particle(numParticles);
    vector<vector<double> > parameters(numParticles);
    vector<vector<int> > exclusionList(numParticles);
    for (int i = 0; i < numParticles; i++) {
        force.getParticleParameters(i, parameters[i]);
        particle[i] = i;
        exclusionList[i].push_back(i);
    }
    for (int i = 0; i < (int)exclusions.size(); i++) {
        exclusionList[exclusions[i].first].push_back(exclusions[i].second);
        exclusionList[exclusions[i].second].push_back(exclusions[i].first);
    }
    CudaNonbondedMethod method = NO_CUTOFF;
    if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
        method = CUTOFF;
    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic) {
        Vec3 boxVectors[3];
        force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
        gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
        method = PERIODIC;
    }
    data.customNonbondedMethod = method;

    // Initialize exceptions.

    int numExceptions = exceptions.size();
    vector<int> exceptionParticle1(numExceptions);
    vector<int> exceptionParticle2(numExceptions);
    vector<vector<double> > exceptionParams(numExceptions);
    for (int i = 0; i < numExceptions; i++)
        force.getExceptionParameters(exceptions[i], exceptionParticle1[i], exceptionParticle2[i], exceptionParams[i]);

436
437
438
439
440
441
442
443
444
445
446
    // Record the tabulated functions.

    for (int i = 0; i < force.getNumFunctions(); i++) {
        string name;
        vector<double> values;
        double min, max;
        bool interpolating;
        force.getFunctionParameters(i, name, values, min, max, interpolating);
        gpuSetTabulatedFunction(gpu, i, name, values, min, max, interpolating);
    }

447
448
449
450
451
452
453
454
    // Record information for the expressions.

    vector<string> paramNames;
    vector<string> combiningRules;
    for (int i = 0; i < force.getNumParameters(); i++) {
        paramNames.push_back(force.getParameterName(i));
        combiningRules.push_back(force.getParameterCombiningRule(i));
    }
455
456
457
458
459
460
    globalParamNames.resize(force.getNumGlobalParameters());
    globalParamValues.resize(force.getNumGlobalParameters());
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
        globalParamNames[i] = force.getGlobalParameterName(i);
        globalParamValues[i] = force.getGlobalParameterDefaultValue(i);
    }
461
    gpuSetCustomNonbondedParameters(gpu, parameters, exclusionList, exceptionParticle1, exceptionParticle2, exceptionParams, method,
462
463
            (float)force.getCutoffDistance(), force.getEnergyFunction(), combiningRules, paramNames, globalParamNames);
    if (globalParamValues.size() > 0)
464
        SetCustomNonbondedGlobalParams(&globalParamValues[0]);
465
466
467
}

void CudaCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
468
469
    if (data.primaryKernel == this) {
        updateGlobalParams(context);
470
        calcForces(context, data);
471
    }
472
473
474
}

double CudaCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
475
476
    if (data.primaryKernel == this) {
        updateGlobalParams(context);
477
        return calcEnergy(context, data, system);
478
    }
479
480
481
    return 0.0;
}

482
483
484
485
486
487
488
489
490
void CudaCalcCustomNonbondedForceKernel::updateGlobalParams(ContextImpl& context) {
    bool changed = false;
    for (int i = 0; i < globalParamNames.size(); i++) {
        float value = (float) context.getParameter(globalParamNames[i]);
        if (value != globalParamValues[i])
            changed = true;
        globalParamValues[i] = value;
    }
    if (changed)
491
        SetCustomNonbondedGlobalParams(&globalParamValues[0]);
492
493
}

494
CudaCalcGBSAOBCForceKernel::~CudaCalcGBSAOBCForceKernel() {
495
496
}

497
void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
498

Peter Eastman's avatar
Peter Eastman committed
499
    int numParticles = system.getNumParticles();
500
    _gpuContext* gpu = data.gpu;
Peter Eastman's avatar
Peter Eastman committed
501
502
    vector<float> radius(numParticles);
    vector<float> scale(numParticles);
503
    vector<float> charge(numParticles);
Peter Eastman's avatar
Peter Eastman committed
504
    for (int i = 0; i < numParticles; i++) {
505
506
        double particleCharge, particleRadius, scalingFactor;
        force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor);
Peter Eastman's avatar
Peter Eastman committed
507
        radius[i] = (float) particleRadius;
508
        scale[i] = (float) scalingFactor;
509
        charge[i] = (float) particleCharge;
510
    }
511
    gpuSetObcParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), radius, scale, charge);
512
513
}

514
void CudaCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
515
516
}

517
static void initializeIntegration(const System& system, CudaPlatform::PlatformData& data, const Integrator& integrator) {
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534

    // Initialize any terms that haven't already been handled by a Force.

    _gpuContext* gpu = data.gpu;
    if (!data.hasBonds)
        gpuSetBondParameters(gpu, vector<int>(), vector<int>(), vector<float>(), vector<float>());
    if (!data.hasAngles)
        gpuSetBondAngleParameters(gpu, vector<int>(), vector<int>(), vector<int>(), vector<float>(), vector<float>());
    if (!data.hasPeriodicTorsions)
        gpuSetDihedralParameters(gpu, vector<int>(), vector<int>(), vector<int>(), vector<int>(), vector<float>(), vector<float>(), vector<int>());
    if (!data.hasRB)
        gpuSetRbDihedralParameters(gpu, vector<int>(), vector<int>(), vector<int>(), vector<int>(), vector<float>(), vector<float>(),
                vector<float>(), vector<float>(), vector<float>(), vector<float>());
    if (!data.hasNonbonded) {
        gpuSetCoulombParameters(gpu, 138.935485f, vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<char>(), vector<vector<int> >(), NO_CUTOFF);
        gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, vector<int>(), vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<float>());
    }
535
536
537
    
    // Set masses.
    
Peter Eastman's avatar
Peter Eastman committed
538
539
540
541
    int numParticles = system.getNumParticles();
    vector<float> mass(numParticles);
    for (int i = 0; i < numParticles; i++)
        mass[i] = (float) system.getParticleMass(i);
542
543
544
545
    gpuSetMass(gpu, mass);
    
    // Set constraints.
    
546
    int numConstraints = system.getNumConstraints();
Peter Eastman's avatar
Peter Eastman committed
547
548
    vector<int> particle1(numConstraints);
    vector<int> particle2(numConstraints);
549
550
551
552
    vector<float> distance(numConstraints);
    vector<float> invMass1(numConstraints);
    vector<float> invMass2(numConstraints);
    for (int i = 0; i < numConstraints; i++) {
Peter Eastman's avatar
Peter Eastman committed
553
        int particle1Index, particle2Index;
554
        double constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
555
556
557
        system.getConstraintParameters(i, particle1Index, particle2Index, constraintDistance);
        particle1[i] = particle1Index;
        particle2[i] = particle2Index;
558
        distance[i] = (float) constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
559
560
        invMass1[i] = 1.0f/mass[particle1Index];
        invMass2[i] = 1.0f/mass[particle2Index];
561
    }
562
    gpuSetConstraintParameters(gpu, particle1, particle2, distance, invMass1, invMass2, (float)integrator.getConstraintTolerance());
563
564
    
    // Finish initialization.
565

566
567
568
    gpuBuildThreadBlockWorkList(gpu);
    gpuBuildExclusionList(gpu);
    gpuBuildOutputBuffers(gpu);
569
    gpuSetConstants(gpu);
570
571
572
    kClearBornForces(gpu);
    kClearForces(gpu);
    cudaThreadSynchronize();
573
574
}

575
double CudaCalcGBSAOBCForceKernel::executeEnergy(ContextImpl& context) {
576
	return 0.0;
577
578
579
580
581
582
}

CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
}

void CudaIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
583
    initializeIntegration(system, data, integrator);
584
585
586
    prevStepSize = -1.0;
}

587
void CudaIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
588
589
590
591
592
    _gpuContext* gpu = data.gpu;
    double stepSize = integrator.getStepSize();
    if (stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
593
        gpuSetVerletIntegrationParameters(gpu, (float) stepSize, 0.0f);
594
595
596
597
598
        gpuSetConstants(gpu);
        prevStepSize = stepSize;
    }
    kVerletUpdatePart1(gpu);
    kApplyFirstShake(gpu);
599
    kApplyFirstSettle(gpu);
600
    kApplyFirstCCMA(gpu);
601
602
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
603
604
            gpu->bCalculateCM = true;
    kVerletUpdatePart2(gpu);
605
    data.time += stepSize;
606
    data.stepCount++;
607
608
609
610
611
612
}

CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
}

void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
613
    initializeIntegration(system, data, integrator);
614
615
616
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
617
618
619
    prevStepSize = -1.0;
}

620
void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
621
    _gpuContext* gpu = data.gpu;
622
623
624
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
625
626
627
628
    if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
        double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
629
        gpuSetLangevinIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature, 0.0f);
630
631
632
633
634
635
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
636
    kLangevinUpdatePart1(gpu);
637
    kApplyFirstShake(gpu);
638
    kApplyFirstSettle(gpu);
639
    kApplyFirstCCMA(gpu);
640
641
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
642
            gpu->bCalculateCM = true;
643
    kLangevinUpdatePart2(gpu);
644
    kApplySecondShake(gpu);
645
    kApplySecondSettle(gpu);
646
    kApplySecondCCMA(gpu);
647
    data.time += stepSize;
648
    data.stepCount++;
649
}
650
651
652
653
654

CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
}

void CudaIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
655
    initializeIntegration(system, data, integrator);
656
657
658
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
659
660
661
    prevStepSize = -1.0;
}

662
void CudaIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
663
664
665
666
667
668
669
670
    _gpuContext* gpu = data.gpu;
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
    if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
        double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
671
        gpuSetBrownianIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature);
672
673
674
675
676
677
678
679
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
    kBrownianUpdatePart1(gpu);
    kApplyFirstShake(gpu);
680
    kApplyFirstSettle(gpu);
681
    kApplyFirstCCMA(gpu);
682
683
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
684
685
            gpu->bCalculateCM = true;
    kBrownianUpdatePart2(gpu);
686
    data.time += stepSize;
687
688
689
690
691
692
693
694
695
696
697
    data.stepCount++;
}

CudaIntegrateVariableVerletStepKernel::~CudaIntegrateVariableVerletStepKernel() {
}

void CudaIntegrateVariableVerletStepKernel::initialize(const System& system, const VariableVerletIntegrator& integrator) {
    initializeIntegration(system, data, integrator);
    prevErrorTol = -1.0;
}

698
void CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
699
700
701
702
703
704
705
706
707
    _gpuContext* gpu = data.gpu;
    double errorTol = integrator.getErrorTolerance();
    if (errorTol != prevErrorTol) {
        // Initialize the GPU parameters.

        gpuSetVerletIntegrationParameters(gpu, 0.0f, (float) errorTol);
        gpuSetConstants(gpu);
        prevErrorTol = errorTol;
    }
708
    float maxStepSize = (float)(maxTime-data.time);
709
    kSelectVerletStepSize(gpu, maxStepSize);
710
711
712
713
714
715
716
717
718
719
    kVerletUpdatePart1(gpu);
    kApplyFirstShake(gpu);
    kApplyFirstSettle(gpu);
    kApplyFirstCCMA(gpu);
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
            gpu->bCalculateCM = true;
    kVerletUpdatePart2(gpu);
    gpu->psStepSize->Download();
    data.time += (*gpu->psStepSize)[0].y;
720
721
    if ((*gpu->psStepSize)[0].y == maxStepSize)
        data.time = maxTime; // Avoid round-off error
722
    data.stepCount++;
723
}
724

725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
CudaIntegrateVariableLangevinStepKernel::~CudaIntegrateVariableLangevinStepKernel() {
}

void CudaIntegrateVariableLangevinStepKernel::initialize(const System& system, const VariableLangevinIntegrator& integrator) {
    initializeIntegration(system, data, integrator);
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
    prevErrorTol = -1.0;
}

void CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
    _gpuContext* gpu = data.gpu;
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double errorTol = integrator.getErrorTolerance();
    if (temperature != prevTemp || friction != prevFriction || errorTol != prevErrorTol) {
        // Initialize the GPU parameters.

        double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
        gpuSetLangevinIntegrationParameters(gpu, (float) tau, 0.0f, (float) temperature, errorTol);
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevErrorTol = errorTol;
    }
    float maxStepSize = (float)(maxTime-data.time);
    kSelectLangevinStepSize(gpu, maxStepSize);
    kLangevinUpdatePart1(gpu);
    kApplyFirstShake(gpu);
    kApplyFirstSettle(gpu);
    kApplyFirstCCMA(gpu);
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
            gpu->bCalculateCM = true;
    kLangevinUpdatePart2(gpu);
    kApplySecondShake(gpu);
    kApplySecondSettle(gpu);
    kApplySecondCCMA(gpu);
    gpu->psStepSize->Download();
    data.time += (*gpu->psStepSize)[0].y;
    if ((*gpu->psStepSize)[0].y == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
}

772
773
774
775
CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
}

void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
776
777
778
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) thermostat.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
779
780
781
    prevStepSize = -1.0;
}

782
void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
783
    _gpuContext* gpu = data.gpu;
784
785
    double temperature = context.getParameter(AndersenThermostat::Temperature());
    double frequency = context.getParameter(AndersenThermostat::CollisionFrequency());
786
787
788
789
    double stepSize = context.getIntegrator().getStepSize();
    if (temperature != prevTemp || frequency != prevFrequency || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
790
        gpuSetAndersenThermostatParameters(gpu, (float) temperature, frequency);
791
792
793
794
795
796
797
798
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFrequency = frequency;
        prevStepSize = stepSize;
    }
    kCalculateAndersenThermostat(gpu);
}
799

800
void CudaCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
801
802
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
803
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
804
        masses[i] = system.getParticleMass(i);
805
806
}

807
double CudaCalcKineticEnergyKernel::execute(ContextImpl& context) {
808
809
810
    // We don't currently have a GPU kernel to do this, so we retrieve the velocities and calculate the energy
    // on the CPU.
    
811
    const Stream& velocities = context.getVelocities();
812
813
814
815
816
817
818
819
    double* v = new double[velocities.getSize()*3];
    velocities.saveToArray(v);
    double energy = 0.0;
    for (size_t i = 0; i < masses.size(); ++i)
        energy += masses[i]*(v[i*3]*v[i*3]+v[i*3+1]*v[i*3+1]+v[i*3+2]*v[i*3+2]);
    delete v;
    return 0.5*energy;
}
820

821
void CudaRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
822
    data.removeCM = true;
823
    data.cmMotionFrequency = force.getFrequency();
824
825
}

826
void CudaRemoveCMMotionKernel::execute(ContextImpl& context) {
827
}