CudaKernels.cpp 31 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, System& system) {
static double calcEnergy(ContextImpl& context, CudaPlatform::PlatformData& data, System& system) {
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
    // New section 2009-09-03: calculate energies and forces, then return reduced energies

    _gpuContext* gpu = data.gpu;
    
    if (gpu->sim.nonbondedMethod == EWALD)
    {
        // We don't currently have GPU kernels to calculate energy, so instead we have the reference
        // platform do it.  This is VERY slow.

        LangevinIntegrator integrator(0.0, 1.0, 0.0);
        ReferencePlatform platform;
        Context refContext(system, integrator, platform);
        const Stream& positions = context.getPositions();
        double* posData = new double[positions.getSize()*3];
        positions.saveToArray(posData);
        vector<Vec3> pos(positions.getSize());
        for (int i = 0; i < (int)pos.size(); i++)
            pos[i] = Vec3(posData[3*i], posData[3*i+1], posData[3*i+2]);
        delete[] posData;
        refContext.setPositions(pos);
 //       printf("Total CPU energies: %f\n", refContext.getState(State::Energy).getPotentialEnergy());
        return refContext.getState(State::Energy).getPotentialEnergy();
    }
    else
    {
        if (data.nonbondedMethod != NO_CUTOFF && data.stepCount%100 == 0)
            gpuReorderAtoms(gpu);
        data.stepCount++;
93
        kClearEnergy(gpu);
94
95
        if (gpu->bIncludeGBSA) {
            gpu->bRecalculateBornRadii = true;
96
97
            kCalculateCDLJObcGbsaForces1(gpu);
            kReduceObcGbsaBornForces(gpu);
98
99
100
            kCalculateObcGbsaForces2(gpu);
        }
        else if (data.hasNonbonded) {
101
            kCalculateCDLJForces(gpu);
102
        }
103
        kCalculateLocalForces(gpu);
104
105
        if (gpu->bIncludeGBSA)
            kReduceBornSumAndForces(gpu);
106
        return kReduceEnergy(gpu);
107
108
    }
    return 0.0f;
109
110
}

111
112
113
void CudaInitializeForcesKernel::initialize(const System& system) {
}

114
void CudaInitializeForcesKernel::execute(ContextImpl& context) {
115
116
}

117
118
119
void CudaUpdateTimeKernel::initialize(const System& system) {
}

120
double CudaUpdateTimeKernel::getTime(const ContextImpl& context) const {
121
122
123
    return data.time;
}

124
void CudaUpdateTimeKernel::setTime(ContextImpl& context, double time) {
125
126
127
    data.time = time;
}

128
129
130
131
132
133
134
CudaCalcHarmonicBondForceKernel::~CudaCalcHarmonicBondForceKernel() {
}

void CudaCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasBonds = true;
135
    numBonds = force.getNumBonds();
Peter Eastman's avatar
Peter Eastman committed
136
137
    vector<int> particle1(numBonds);
    vector<int> particle2(numBonds);
138
139
140
141
    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
142
        force.getBondParameters(i, particle1[i], particle2[i], lengthValue, kValue);
143
144
145
        length[i] = (float) lengthValue;
        k[i] = (float) kValue;
    }
Peter Eastman's avatar
Peter Eastman committed
146
    gpuSetBondParameters(data.gpu, particle1, particle2, length, k);
147
148
}

149
void CudaCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
150
151
152
153
    if (data.primaryKernel == this)
        calcForces(context, data);
}

154
double CudaCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
155
    if (data.primaryKernel == this)
156
        return calcEnergy(context, data, system);
157
158
159
160
161
162
163
164
165
166
    return 0.0;
}

CudaCalcHarmonicAngleForceKernel::~CudaCalcHarmonicAngleForceKernel() {
}

void CudaCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasAngles = true;
167
    numAngles = force.getNumAngles();
168
    const float RadiansToDegrees = (float) (180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
169
170
171
    vector<int> particle1(numAngles);
    vector<int> particle2(numAngles);
    vector<int> particle3(numAngles);
172
173
174
175
    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
176
        force.getAngleParameters(i, particle1[i], particle2[i], particle3[i], angleValue, kValue);
177
178
        angle[i] = (float) (angleValue*RadiansToDegrees);
        k[i] = (float) kValue;
179
    }
Peter Eastman's avatar
Peter Eastman committed
180
    gpuSetBondAngleParameters(data.gpu, particle1, particle2, particle3, angle, k);
181
}
182

183
void CudaCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
184
185
186
187
    if (data.primaryKernel == this)
        calcForces(context, data);
}

188
double CudaCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
189
    if (data.primaryKernel == this)
190
        return calcEnergy(context, data, system);
191
192
193
194
195
196
197
198
199
200
201
    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();
202
    const float RadiansToDegrees = (float)(180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
203
204
205
206
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
207
208
209
210
211
    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
212
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], periodicity[i], phaseValue, kValue);
213
214
        k[i] = (float) kValue;
        phase[i] = (float) (phaseValue*RadiansToDegrees);
215
    }
Peter Eastman's avatar
Peter Eastman committed
216
    gpuSetDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, k, phase, periodicity);
217
218
}

219
void CudaCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
220
221
222
223
    if (data.primaryKernel == this)
        calcForces(context, data);
}

224
double CudaCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
225
    if (data.primaryKernel == this)
226
        return calcEnergy(context, data, system);
227
228
229
230
231
232
233
234
235
236
237
    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
238
239
240
241
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
242
243
244
245
246
247
248
249
    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
250
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], c[0], c[1], c[2], c[3], c[4], c[5]);
251
252
253
254
255
256
        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];
257
    }
Peter Eastman's avatar
Peter Eastman committed
258
    gpuSetRbDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
259
260
}

261
void CudaCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
262
263
264
265
    if (data.primaryKernel == this)
        calcForces(context, data);
}

266
double CudaCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
267
    if (data.primaryKernel == this)
268
        return calcEnergy(context, data, system);
269
270
271
272
273
274
    return 0.0;
}

CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() {
}

275
void CudaCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
276
277
278
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasNonbonded = true;
Peter Eastman's avatar
Peter Eastman committed
279
    numParticles = force.getNumParticles();
280
    _gpuContext* gpu = data.gpu;
281
282
283
284

    // Identify which exceptions are 1-4 interactions.

    vector<pair<int, int> > exclusions;
285
    vector<int> exceptions;
286
287
288
289
290
291
    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)
292
            exceptions.push_back(i);
293
294
    }

295
296
    // Initialize nonbonded interactions.
    
297
    {
Peter Eastman's avatar
Peter Eastman committed
298
299
300
301
        vector<int> particle(numParticles);
        vector<float> c6(numParticles);
        vector<float> c12(numParticles);
        vector<float> q(numParticles);
302
        vector<char> symbol;
Peter Eastman's avatar
Peter Eastman committed
303
304
        vector<vector<int> > exclusionList(numParticles);
        for (int i = 0; i < numParticles; i++) {
305
            double charge, radius, depth;
Peter Eastman's avatar
Peter Eastman committed
306
307
            force.getParticleParameters(i, charge, radius, depth);
            particle[i] = i;
308
309
310
            q[i] = (float) charge;
            c6[i] = (float) (4*depth*pow(radius, 6.0));
            c12[i] = (float) (4*depth*pow(radius, 12.0));
311
312
            exclusionList[i].push_back(i);
        }
313
        for (int i = 0; i < (int)exclusions.size(); i++) {
314
315
316
            exclusionList[exclusions[i].first].push_back(exclusions[i].second);
            exclusionList[exclusions[i].second].push_back(exclusions[i].first);
        }
317
318
        CudaNonbondedMethod method = NO_CUTOFF;
        if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
319
            gpuSetNonbondedCutoff(gpu, (float)force.getCutoffDistance(), force.getReactionFieldDielectric());
320
321
322
323
324
            method = CUTOFF;
        }
        if (force.getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
            Vec3 boxVectors[3];
            force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
325
            gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
326
327
            method = PERIODIC;
        }
328
329
330
331

        if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
            Vec3 boxVectors[3];
            force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
332
            gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
333
334
335
336
337
338
            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;
339
340
341
            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));
342
343
344
345
346
347
            if (kmaxx%2 == 0)
                kmaxx++;
            if (kmaxy%2 == 0)
                kmaxy++;
            if (kmaxz%2 == 0)
                kmaxz++;
348
            gpuSetEwaldParameters(gpu, (float)alpha, kmaxx, kmaxy, kmaxz);
349
350
            method = EWALD;
        }
351
352
        data.nonbondedMethod = method;
        gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method);
353
354
355
356
    }

    // Initialize 1-4 nonbonded interactions.
    
357
    {
358
359
360
361
362
363
364
365
        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++) {
366
            double charge, sig, eps;
367
            force.getExceptionParameters(exceptions[i], particle1[i], particle2[i], charge, sig, eps);
368
369
            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
370
371
            q1[i] = (float) charge;
            q2[i] = 1.0f;
372
        }
Peter Eastman's avatar
Peter Eastman committed
373
        gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, particle1, particle2, c6, c12, q1, q2);
374
375
376
    }
}

377
void CudaCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
378
379
    if (data.primaryKernel == this)
        calcForces(context, data);
380
381
}

382
double CudaCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
383
    if (data.primaryKernel == this)
384
        return calcEnergy(context, data, system);
385
    return 0.0;
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
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
CudaCalcCustomNonbondedForceKernel::~CudaCalcCustomNonbondedForceKernel() {
}

void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    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]);

    // 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));
    }
    gpuSetCustomNonbondedParameters(gpu, parameters, exclusionList, exceptionParticle1, exceptionParticle2, exceptionParams, method,
            (float)force.getCutoffDistance(), force.getEnergyFunction(), combiningRules, paramNames);
}

void CudaCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
    if (data.primaryKernel == this)
        calcForces(context, data);
}

double CudaCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
    if (data.primaryKernel == this)
466
        return calcEnergy(context, data, system);
467
468
469
    return 0.0;
}

470
CudaCalcGBSAOBCForceKernel::~CudaCalcGBSAOBCForceKernel() {
471
472
}

473
void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
474

Peter Eastman's avatar
Peter Eastman committed
475
    int numParticles = system.getNumParticles();
476
    _gpuContext* gpu = data.gpu;
Peter Eastman's avatar
Peter Eastman committed
477
478
    vector<float> radius(numParticles);
    vector<float> scale(numParticles);
479
    vector<float> charge(numParticles);
Peter Eastman's avatar
Peter Eastman committed
480
    for (int i = 0; i < numParticles; i++) {
481
482
        double particleCharge, particleRadius, scalingFactor;
        force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor);
Peter Eastman's avatar
Peter Eastman committed
483
        radius[i] = (float) particleRadius;
484
        scale[i] = (float) scalingFactor;
485
        charge[i] = (float) particleCharge;
486
    }
487
    gpuSetObcParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), radius, scale, charge);
488
489
}

490
void CudaCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
491
492
}

493
static void initializeIntegration(const System& system, CudaPlatform::PlatformData& data, const Integrator& integrator) {
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510

    // 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>());
    }
511
512
513
    
    // Set masses.
    
Peter Eastman's avatar
Peter Eastman committed
514
515
516
517
    int numParticles = system.getNumParticles();
    vector<float> mass(numParticles);
    for (int i = 0; i < numParticles; i++)
        mass[i] = (float) system.getParticleMass(i);
518
519
520
521
    gpuSetMass(gpu, mass);
    
    // Set constraints.
    
522
    int numConstraints = system.getNumConstraints();
Peter Eastman's avatar
Peter Eastman committed
523
524
    vector<int> particle1(numConstraints);
    vector<int> particle2(numConstraints);
525
526
527
528
    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
529
        int particle1Index, particle2Index;
530
        double constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
531
532
533
        system.getConstraintParameters(i, particle1Index, particle2Index, constraintDistance);
        particle1[i] = particle1Index;
        particle2[i] = particle2Index;
534
        distance[i] = (float) constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
535
536
        invMass1[i] = 1.0f/mass[particle1Index];
        invMass2[i] = 1.0f/mass[particle2Index];
537
    }
538
    gpuSetConstraintParameters(gpu, particle1, particle2, distance, invMass1, invMass2, (float)integrator.getConstraintTolerance());
539
540
    
    // Finish initialization.
541

542
543
544
    gpuBuildThreadBlockWorkList(gpu);
    gpuBuildExclusionList(gpu);
    gpuBuildOutputBuffers(gpu);
545
    gpuSetConstants(gpu);
546
547
548
    kClearBornForces(gpu);
    kClearForces(gpu);
    cudaThreadSynchronize();
549
550
}

551
double CudaCalcGBSAOBCForceKernel::executeEnergy(ContextImpl& context) {
552
	return 0.0;
553
554
555
556
557
558
}

CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
}

void CudaIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
559
    initializeIntegration(system, data, integrator);
560
561
562
    prevStepSize = -1.0;
}

563
void CudaIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
564
565
566
567
568
    _gpuContext* gpu = data.gpu;
    double stepSize = integrator.getStepSize();
    if (stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
569
        gpuSetVerletIntegrationParameters(gpu, (float) stepSize, 0.0f);
570
571
572
573
574
        gpuSetConstants(gpu);
        prevStepSize = stepSize;
    }
    kVerletUpdatePart1(gpu);
    kApplyFirstShake(gpu);
575
    kApplyFirstSettle(gpu);
576
    kApplyFirstCCMA(gpu);
577
578
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
579
580
            gpu->bCalculateCM = true;
    kVerletUpdatePart2(gpu);
581
    data.time += stepSize;
582
    data.stepCount++;
583
584
585
586
587
588
}

CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
}

void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
589
    initializeIntegration(system, data, integrator);
590
591
592
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
593
594
595
    prevStepSize = -1.0;
}

596
void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
597
    _gpuContext* gpu = data.gpu;
598
599
600
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
601
602
603
604
    if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
        double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
605
        gpuSetLangevinIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature, 0.0f);
606
607
608
609
610
611
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
612
    kLangevinUpdatePart1(gpu);
613
    kApplyFirstShake(gpu);
614
    kApplyFirstSettle(gpu);
615
    kApplyFirstCCMA(gpu);
616
617
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
618
            gpu->bCalculateCM = true;
619
    kLangevinUpdatePart2(gpu);
620
    kApplySecondShake(gpu);
621
    kApplySecondSettle(gpu);
622
    kApplySecondCCMA(gpu);
623
    data.time += stepSize;
624
    data.stepCount++;
625
}
626
627
628
629
630

CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
}

void CudaIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
631
    initializeIntegration(system, data, integrator);
632
633
634
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
635
636
637
    prevStepSize = -1.0;
}

638
void CudaIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
639
640
641
642
643
644
645
646
    _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);
647
        gpuSetBrownianIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature);
648
649
650
651
652
653
654
655
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
    kBrownianUpdatePart1(gpu);
    kApplyFirstShake(gpu);
656
    kApplyFirstSettle(gpu);
657
    kApplyFirstCCMA(gpu);
658
659
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
660
661
            gpu->bCalculateCM = true;
    kBrownianUpdatePart2(gpu);
662
    data.time += stepSize;
663
664
665
666
667
668
669
670
671
672
673
    data.stepCount++;
}

CudaIntegrateVariableVerletStepKernel::~CudaIntegrateVariableVerletStepKernel() {
}

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

674
void CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
675
676
677
678
679
680
681
682
683
    _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;
    }
684
    float maxStepSize = (float)(maxTime-data.time);
685
    kSelectVerletStepSize(gpu, maxStepSize);
686
687
688
689
690
691
692
693
694
695
    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;
696
697
    if ((*gpu->psStepSize)[0].y == maxStepSize)
        data.time = maxTime; // Avoid round-off error
698
    data.stepCount++;
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
736
737
738
739
740
741
742
743
744
745
746
747
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++;
}

748
749
750
751
CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
}

void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
752
753
754
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) thermostat.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
755
756
757
    prevStepSize = -1.0;
}

758
void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
759
    _gpuContext* gpu = data.gpu;
760
761
    double temperature = context.getParameter(AndersenThermostat::Temperature());
    double frequency = context.getParameter(AndersenThermostat::CollisionFrequency());
762
763
764
765
    double stepSize = context.getIntegrator().getStepSize();
    if (temperature != prevTemp || frequency != prevFrequency || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
766
        gpuSetAndersenThermostatParameters(gpu, (float) temperature, frequency);
767
768
769
770
771
772
773
774
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFrequency = frequency;
        prevStepSize = stepSize;
    }
    kCalculateAndersenThermostat(gpu);
}
775

776
void CudaCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
777
778
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
779
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
780
        masses[i] = system.getParticleMass(i);
781
782
}

783
double CudaCalcKineticEnergyKernel::execute(ContextImpl& context) {
784
785
786
    // We don't currently have a GPU kernel to do this, so we retrieve the velocities and calculate the energy
    // on the CPU.
    
787
    const Stream& velocities = context.getVelocities();
788
789
790
791
792
793
794
795
    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;
}
796

797
void CudaRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
798
    data.removeCM = true;
799
    data.cmMotionFrequency = force.getFrequency();
800
801
}

802
void CudaRemoveCMMotionKernel::execute(ContextImpl& context) {
803
}