CudaKernels.cpp 30.7 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
31
#include "openmm/Context.h"
#include "openmm/internal/ContextImpl.h"
32
#include "kernels/gputypes.h"
33
#include "kernels/cudaKernels.h"
34
35
36
37
38
39
40
#include <cmath>

extern "C" int gpuSetConstants( gpuContext gpu );

using namespace OpenMM;
using namespace std;

41
42
43
44
void CudaCalcForcesAndEnergyKernel::initialize(const System& system) {
}

void CudaCalcForcesAndEnergyKernel::beginForceComputation(ContextImpl& context) {
45
    _gpuContext* gpu = data.gpu;
46
    if (data.nonbondedMethod != NO_CUTOFF && data.computeForceCount%100 == 0)
47
        gpuReorderAtoms(gpu);
48
    data.computeForceCount++;
49
    kClearForces(gpu);
50
51
52
53
}

void CudaCalcForcesAndEnergyKernel::finishForceComputation(ContextImpl& context) {
    _gpuContext* gpu = data.gpu;
54
    if (gpu->bIncludeGBSA) {
55
        gpu->bRecalculateBornRadii = true;
56
57
58
59
        kCalculateCDLJObcGbsaForces1(gpu);
        kReduceObcGbsaBornForces(gpu);
        kCalculateObcGbsaForces2(gpu);
    }
60
    else if (data.hasNonbonded)
61
        kCalculateCDLJForces(gpu);
62
63
    if (data.hasCustomNonbonded)
        kCalculateCustomNonbondedForces(gpu, data.hasNonbonded);
64
    kCalculateLocalForces(gpu);
65
    kReduceForces(gpu);
66
67
}

68
void CudaCalcForcesAndEnergyKernel::beginEnergyComputation(ContextImpl& context) {
69
    _gpuContext* gpu = data.gpu;
70
71
72
73
    if (data.nonbondedMethod != NO_CUTOFF && data.stepCount%100 == 0)
        gpuReorderAtoms(gpu);
    data.stepCount++;
    kClearEnergy(gpu);
74
75
76
77
}

double CudaCalcForcesAndEnergyKernel::finishEnergyComputation(ContextImpl& context) {
    _gpuContext* gpu = data.gpu;
78
79
80
81
82
    if (gpu->bIncludeGBSA) {
        gpu->bRecalculateBornRadii = true;
        kCalculateCDLJObcGbsaForces1(gpu);
        kReduceObcGbsaBornForces(gpu);
        kCalculateObcGbsaForces2(gpu);
83
    }
84
85
86
87
88
89
90
91
    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;
92
93
}

94
95
96
void CudaUpdateTimeKernel::initialize(const System& system) {
}

97
double CudaUpdateTimeKernel::getTime(const ContextImpl& context) const {
98
99
100
    return data.time;
}

101
void CudaUpdateTimeKernel::setTime(ContextImpl& context, double time) {
102
103
104
    data.time = time;
}

105
106
107
108
109
CudaCalcHarmonicBondForceKernel::~CudaCalcHarmonicBondForceKernel() {
}

void CudaCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
    data.hasBonds = true;
110
    numBonds = force.getNumBonds();
Peter Eastman's avatar
Peter Eastman committed
111
112
    vector<int> particle1(numBonds);
    vector<int> particle2(numBonds);
113
114
115
116
    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
117
        force.getBondParameters(i, particle1[i], particle2[i], lengthValue, kValue);
118
119
120
        length[i] = (float) lengthValue;
        k[i] = (float) kValue;
    }
Peter Eastman's avatar
Peter Eastman committed
121
    gpuSetBondParameters(data.gpu, particle1, particle2, length, k);
122
123
}

124
void CudaCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
125
126
}

127
double CudaCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
128
129
130
131
132
133
134
135
    return 0.0;
}

CudaCalcHarmonicAngleForceKernel::~CudaCalcHarmonicAngleForceKernel() {
}

void CudaCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
    data.hasAngles = true;
136
    numAngles = force.getNumAngles();
137
    const float RadiansToDegrees = (float) (180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
138
139
140
    vector<int> particle1(numAngles);
    vector<int> particle2(numAngles);
    vector<int> particle3(numAngles);
141
142
143
144
    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
145
        force.getAngleParameters(i, particle1[i], particle2[i], particle3[i], angleValue, kValue);
146
147
        angle[i] = (float) (angleValue*RadiansToDegrees);
        k[i] = (float) kValue;
148
    }
Peter Eastman's avatar
Peter Eastman committed
149
    gpuSetBondAngleParameters(data.gpu, particle1, particle2, particle3, angle, k);
150
}
151

152
void CudaCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
153
154
}

155
double CudaCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
156
157
158
159
160
161
162
163
164
    return 0.0;
}

CudaCalcPeriodicTorsionForceKernel::~CudaCalcPeriodicTorsionForceKernel() {
}

void CudaCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
    data.hasPeriodicTorsions = true;
    numTorsions = force.getNumTorsions();
165
    const float RadiansToDegrees = (float)(180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
166
167
168
169
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
170
171
172
173
174
    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
175
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], periodicity[i], phaseValue, kValue);
176
177
        k[i] = (float) kValue;
        phase[i] = (float) (phaseValue*RadiansToDegrees);
178
    }
Peter Eastman's avatar
Peter Eastman committed
179
    gpuSetDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, k, phase, periodicity);
180
181
}

182
void CudaCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
183
184
}

185
double CudaCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
186
187
188
189
190
191
192
193
194
    return 0.0;
}

CudaCalcRBTorsionForceKernel::~CudaCalcRBTorsionForceKernel() {
}

void CudaCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
    data.hasRB = true;
    numTorsions = force.getNumTorsions();
Peter Eastman's avatar
Peter Eastman committed
195
196
197
198
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
199
200
201
202
203
204
205
206
    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
207
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], c[0], c[1], c[2], c[3], c[4], c[5]);
208
209
210
211
212
213
        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];
214
    }
Peter Eastman's avatar
Peter Eastman committed
215
    gpuSetRbDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
216
217
}

218
void CudaCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
219
220
}

221
double CudaCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
222
223
224
225
226
227
    return 0.0;
}

CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() {
}

228
void CudaCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
229
    data.hasNonbonded = true;
Peter Eastman's avatar
Peter Eastman committed
230
    numParticles = force.getNumParticles();
231
    _gpuContext* gpu = data.gpu;
232
233
234
235

    // Identify which exceptions are 1-4 interactions.

    vector<pair<int, int> > exclusions;
236
    vector<int> exceptions;
237
238
239
240
241
242
    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)
243
            exceptions.push_back(i);
244
245
    }

246
247
    // Initialize nonbonded interactions.
    
248
    {
Peter Eastman's avatar
Peter Eastman committed
249
250
251
252
        vector<int> particle(numParticles);
        vector<float> c6(numParticles);
        vector<float> c12(numParticles);
        vector<float> q(numParticles);
253
        vector<char> symbol;
Peter Eastman's avatar
Peter Eastman committed
254
255
        vector<vector<int> > exclusionList(numParticles);
        for (int i = 0; i < numParticles; i++) {
256
            double charge, radius, depth;
Peter Eastman's avatar
Peter Eastman committed
257
258
            force.getParticleParameters(i, charge, radius, depth);
            particle[i] = i;
259
260
261
            q[i] = (float) charge;
            c6[i] = (float) (4*depth*pow(radius, 6.0));
            c12[i] = (float) (4*depth*pow(radius, 12.0));
262
263
            exclusionList[i].push_back(i);
        }
264
        for (int i = 0; i < (int)exclusions.size(); i++) {
265
266
267
            exclusionList[exclusions[i].first].push_back(exclusions[i].second);
            exclusionList[exclusions[i].second].push_back(exclusions[i].first);
        }
Peter Eastman's avatar
Peter Eastman committed
268
        Vec3 boxVectors[3];
269
        system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
Peter Eastman's avatar
Peter Eastman committed
270
        gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
271
272
        CudaNonbondedMethod method = NO_CUTOFF;
        if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
273
            gpuSetNonbondedCutoff(gpu, (float)force.getCutoffDistance(), force.getReactionFieldDielectric());
274
275
276
277
278
            method = CUTOFF;
        }
        if (force.getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
            method = PERIODIC;
        }
279
        if (force.getNonbondedMethod() == NonbondedForce::Ewald || force.getNonbondedMethod() == NonbondedForce::PME) {
280
281
282
283
284
285
            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
286
287
288
            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));
289
290
291
292
293
294
295
296
297
298
299
            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 {
300
301
302
                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
303
                gpuSetPMEParameters(gpu, (float) alpha, gridSizeX, gridSizeY, gridSizeZ);
304
305
                method = PARTICLE_MESH_EWALD;
            }
306
        }
307
308
        data.nonbondedMethod = method;
        gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method);
309
310
311
312

        // Compute the Ewald self energy.

        data.ewaldSelfEnergy = 0.0;
313
314
315
316
317
        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];
        }
318
319
320
321
    }

    // Initialize 1-4 nonbonded interactions.
    
322
    {
323
324
325
326
327
328
329
330
        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++) {
331
            double charge, sig, eps;
332
            force.getExceptionParameters(exceptions[i], particle1[i], particle2[i], charge, sig, eps);
333
334
            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
335
336
            q1[i] = (float) charge;
            q2[i] = 1.0f;
337
        }
Peter Eastman's avatar
Peter Eastman committed
338
        gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, particle1, particle2, c6, c12, q1, q2);
339
340
341
    }
}

342
void CudaCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
343
344
}

345
double CudaCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
346
    return 0.0;
347
348
}

349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
CudaCalcCustomNonbondedForceKernel::~CudaCalcCustomNonbondedForceKernel() {
}

void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
    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);
    }
386
387
388
    Vec3 boxVectors[3];
    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
    gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
    CudaNonbondedMethod method = NO_CUTOFF;
    if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
        method = CUTOFF;
    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic) {
        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]);

406
407
408
409
410
411
412
413
414
415
416
    // 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);
    }

417
418
419
420
421
422
423
424
    // 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));
    }
425
426
427
428
429
430
    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);
    }
431
    gpuSetCustomNonbondedParameters(gpu, parameters, exclusionList, exceptionParticle1, exceptionParticle2, exceptionParams, method,
432
433
            (float)force.getCutoffDistance(), force.getEnergyFunction(), combiningRules, paramNames, globalParamNames);
    if (globalParamValues.size() > 0)
434
        SetCustomNonbondedGlobalParams(&globalParamValues[0]);
435
436
437
}

void CudaCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
438
    updateGlobalParams(context);
439
440
441
}

double CudaCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
442
    updateGlobalParams(context);
443
444
445
    return 0.0;
}

446
447
448
449
450
451
452
453
454
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)
455
        SetCustomNonbondedGlobalParams(&globalParamValues[0]);
456
457
}

458
CudaCalcGBSAOBCForceKernel::~CudaCalcGBSAOBCForceKernel() {
459
460
}

461
void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
462

Peter Eastman's avatar
Peter Eastman committed
463
    int numParticles = system.getNumParticles();
464
    _gpuContext* gpu = data.gpu;
Peter Eastman's avatar
Peter Eastman committed
465
466
    vector<float> radius(numParticles);
    vector<float> scale(numParticles);
467
    vector<float> charge(numParticles);
Peter Eastman's avatar
Peter Eastman committed
468
    for (int i = 0; i < numParticles; i++) {
469
470
        double particleCharge, particleRadius, scalingFactor;
        force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor);
Peter Eastman's avatar
Peter Eastman committed
471
        radius[i] = (float) particleRadius;
472
        scale[i] = (float) scalingFactor;
473
        charge[i] = (float) particleCharge;
474
    }
475
    gpuSetObcParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), radius, scale, charge);
476
477
}

478
void CudaCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
479
480
}

481
static void initializeIntegration(const System& system, CudaPlatform::PlatformData& data, const Integrator& integrator) {
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498

    // 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>());
    }
499
500
501
    
    // Set masses.
    
Peter Eastman's avatar
Peter Eastman committed
502
503
504
505
    int numParticles = system.getNumParticles();
    vector<float> mass(numParticles);
    for (int i = 0; i < numParticles; i++)
        mass[i] = (float) system.getParticleMass(i);
506
507
508
509
    gpuSetMass(gpu, mass);
    
    // Set constraints.
    
510
    int numConstraints = system.getNumConstraints();
Peter Eastman's avatar
Peter Eastman committed
511
512
    vector<int> particle1(numConstraints);
    vector<int> particle2(numConstraints);
513
514
515
516
    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
517
        int particle1Index, particle2Index;
518
        double constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
519
520
521
        system.getConstraintParameters(i, particle1Index, particle2Index, constraintDistance);
        particle1[i] = particle1Index;
        particle2[i] = particle2Index;
522
        distance[i] = (float) constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
523
524
        invMass1[i] = 1.0f/mass[particle1Index];
        invMass2[i] = 1.0f/mass[particle2Index];
525
    }
526
    gpuSetConstraintParameters(gpu, particle1, particle2, distance, invMass1, invMass2, (float)integrator.getConstraintTolerance());
527
528
    
    // Finish initialization.
529

530
531
532
    gpuBuildThreadBlockWorkList(gpu);
    gpuBuildExclusionList(gpu);
    gpuBuildOutputBuffers(gpu);
533
    gpuSetConstants(gpu);
534
535
536
    kClearBornForces(gpu);
    kClearForces(gpu);
    cudaThreadSynchronize();
537
538
}

539
double CudaCalcGBSAOBCForceKernel::executeEnergy(ContextImpl& context) {
540
	return 0.0;
541
542
543
544
545
546
}

CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
}

void CudaIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
547
    initializeIntegration(system, data, integrator);
548
549
550
    prevStepSize = -1.0;
}

551
void CudaIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
552
553
554
555
556
    _gpuContext* gpu = data.gpu;
    double stepSize = integrator.getStepSize();
    if (stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
557
        gpuSetVerletIntegrationParameters(gpu, (float) stepSize, 0.0f);
558
559
560
561
562
        gpuSetConstants(gpu);
        prevStepSize = stepSize;
    }
    kVerletUpdatePart1(gpu);
    kApplyFirstShake(gpu);
563
    kApplyFirstSettle(gpu);
564
    kApplyFirstCCMA(gpu);
565
566
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
567
568
            gpu->bCalculateCM = true;
    kVerletUpdatePart2(gpu);
569
    data.time += stepSize;
570
    data.stepCount++;
571
572
573
574
575
576
}

CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
}

void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
577
    initializeIntegration(system, data, integrator);
578
579
580
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
581
582
583
    prevStepSize = -1.0;
}

584
void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
585
    _gpuContext* gpu = data.gpu;
586
587
588
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
589
590
591
592
    if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
        double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
593
        gpuSetLangevinIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature, 0.0f);
594
595
596
597
598
599
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
600
    kLangevinUpdatePart1(gpu);
601
    kApplyFirstShake(gpu);
602
    kApplyFirstSettle(gpu);
603
    kApplyFirstCCMA(gpu);
604
605
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
606
            gpu->bCalculateCM = true;
607
    kLangevinUpdatePart2(gpu);
608
    kApplySecondShake(gpu);
609
    kApplySecondSettle(gpu);
610
    kApplySecondCCMA(gpu);
611
    data.time += stepSize;
612
    data.stepCount++;
613
}
614
615
616
617
618

CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
}

void CudaIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
619
    initializeIntegration(system, data, integrator);
620
621
622
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
623
624
625
    prevStepSize = -1.0;
}

626
void CudaIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
627
628
629
630
631
632
633
634
    _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);
635
        gpuSetBrownianIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature);
636
637
638
639
640
641
642
643
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
    kBrownianUpdatePart1(gpu);
    kApplyFirstShake(gpu);
644
    kApplyFirstSettle(gpu);
645
    kApplyFirstCCMA(gpu);
646
647
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
648
649
            gpu->bCalculateCM = true;
    kBrownianUpdatePart2(gpu);
650
    data.time += stepSize;
651
652
653
654
655
656
657
658
659
660
661
    data.stepCount++;
}

CudaIntegrateVariableVerletStepKernel::~CudaIntegrateVariableVerletStepKernel() {
}

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

662
void CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
663
664
665
666
667
668
669
670
671
    _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;
    }
672
    float maxStepSize = (float)(maxTime-data.time);
673
    kSelectVerletStepSize(gpu, maxStepSize);
674
675
676
677
678
679
680
681
682
683
    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;
684
685
    if ((*gpu->psStepSize)[0].y == maxStepSize)
        data.time = maxTime; // Avoid round-off error
686
    data.stepCount++;
687
}
688

689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
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++;
}

736
737
738
739
CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
}

void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
740
741
742
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) thermostat.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
743
744
745
    prevStepSize = -1.0;
}

746
void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
747
    _gpuContext* gpu = data.gpu;
748
749
    double temperature = context.getParameter(AndersenThermostat::Temperature());
    double frequency = context.getParameter(AndersenThermostat::CollisionFrequency());
750
751
752
753
    double stepSize = context.getIntegrator().getStepSize();
    if (temperature != prevTemp || frequency != prevFrequency || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
754
        gpuSetAndersenThermostatParameters(gpu, (float) temperature, frequency);
755
756
757
758
759
760
761
762
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFrequency = frequency;
        prevStepSize = stepSize;
    }
    kCalculateAndersenThermostat(gpu);
}
763

764
void CudaCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
765
766
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
767
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
768
        masses[i] = system.getParticleMass(i);
769
770
}

771
double CudaCalcKineticEnergyKernel::execute(ContextImpl& context) {
772
773
774
    // We don't currently have a GPU kernel to do this, so we retrieve the velocities and calculate the energy
    // on the CPU.
    
775
    const Stream& velocities = context.getVelocities();
776
777
778
779
780
781
782
783
    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;
}
784

785
void CudaRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
786
    data.removeCM = true;
787
    data.cmMotionFrequency = force.getFrequency();
788
789
}

790
void CudaRemoveCMMotionKernel::execute(ContextImpl& context) {
791
}