CudaKernels.cpp 22.9 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
30
#include "openmm/LangevinIntegrator.h"
#include "openmm/OpenMMContext.h"
31
#include "ReferencePlatform.h"
32
#include "openmm/internal/OpenMMContextImpl.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
43
static void calcForces(OpenMMContextImpl& context, CudaPlatform::PlatformData& data) {
    _gpuContext* gpu = data.gpu;
44
    if (data.nonbondedMethod != NO_CUTOFF && data.stepCount%100 == 0)
45
46
47
48
        gpuReorderAtoms(gpu);
    data.stepCount++;
    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
56
57
        kCalculateCDLJForces(gpu);
    }
    kCalculateLocalForces(gpu);
58
/*
59
60
61
    if (gpu->bIncludeGBSA)
        kReduceBornSumAndForces(gpu);
    else
62
*/
63
        kReduceForces(gpu);
64
65
}

66
67
68
static double calcEnergy(OpenMMContextImpl& context, System& system) {
    // We don't currently have GPU kernels to calculate energy, so instead we have the reference
    // platform do it.  This is VERY slow.
69

70
71
72
73
74
75
76
    LangevinIntegrator integrator(0.0, 1.0, 0.0);
    ReferencePlatform platform;
    OpenMMContext refContext(system, integrator, platform);
    const Stream& positions = context.getPositions();
    double* posData = new double[positions.getSize()*3];
    positions.saveToArray(posData);
    vector<Vec3> pos(positions.getSize());
77
    for (int i = 0; i < (int)pos.size(); i++)
78
79
80
81
82
83
        pos[i] = Vec3(posData[3*i], posData[3*i+1], posData[3*i+2]);
    delete[] posData;
    refContext.setPositions(pos);
    return refContext.getState(State::Energy).getPotentialEnergy();
}

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

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

90
91
92
93
94
95
96
97
98
99
100
void CudaUpdateTimeKernel::initialize(const System& system) {
}

double CudaUpdateTimeKernel::getTime(const OpenMMContextImpl& context) const {
    return data.time;
}

void CudaUpdateTimeKernel::setTime(OpenMMContextImpl& context, double time) {
    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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
}

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

double CudaCalcHarmonicBondForceKernel::executeEnergy(OpenMMContextImpl& context) {
    if (data.primaryKernel == this)
        return calcEnergy(context, system);
    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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
void CudaCalcHarmonicAngleForceKernel::executeForces(OpenMMContextImpl& context) {
    if (data.primaryKernel == this)
        calcForces(context, data);
}

double CudaCalcHarmonicAngleForceKernel::executeEnergy(OpenMMContextImpl& context) {
    if (data.primaryKernel == this)
        return calcEnergy(context, system);
    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
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
}

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

double CudaCalcPeriodicTorsionForceKernel::executeEnergy(OpenMMContextImpl& context) {
    if (data.primaryKernel == this)
        return calcEnergy(context, system);
    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
235
236
237
238
239
240
241
242
243
244
245
246
247
}

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

double CudaCalcRBTorsionForceKernel::executeEnergy(OpenMMContextImpl& context) {
    if (data.primaryKernel == this)
        return calcEnergy(context, system);
    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
258
259
260
261
262
263
264
265
266
267
268

    // Identify which exceptions are 1-4 interactions.

    vector<pair<int, int> > exclusions;
    vector<int> nb14s;
    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)
            nb14s.push_back(i);
    }
    num14 = nb14s.size();

269
270
    // Initialize nonbonded interactions.
    
271
    {
Peter Eastman's avatar
Peter Eastman committed
272
273
274
275
        vector<int> particle(numParticles);
        vector<float> c6(numParticles);
        vector<float> c12(numParticles);
        vector<float> q(numParticles);
276
        vector<char> symbol;
Peter Eastman's avatar
Peter Eastman committed
277
278
        vector<vector<int> > exclusionList(numParticles);
        for (int i = 0; i < numParticles; i++) {
279
            double charge, radius, depth;
Peter Eastman's avatar
Peter Eastman committed
280
281
            force.getParticleParameters(i, charge, radius, depth);
            particle[i] = i;
282
283
284
            q[i] = (float) charge;
            c6[i] = (float) (4*depth*pow(radius, 6.0));
            c12[i] = (float) (4*depth*pow(radius, 12.0));
285
286
            exclusionList[i].push_back(i);
        }
287
        for (int i = 0; i < (int)exclusions.size(); i++) {
288
289
290
            exclusionList[exclusions[i].first].push_back(exclusions[i].second);
            exclusionList[exclusions[i].second].push_back(exclusions[i].first);
        }
291
292
        CudaNonbondedMethod method = NO_CUTOFF;
        if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
293
            gpuSetNonbondedCutoff(gpu, (float)force.getCutoffDistance(), 78.3f);
294
295
296
297
298
            method = CUTOFF;
        }
        if (force.getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
            Vec3 boxVectors[3];
            force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
299
            gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
300
301
            method = PERIODIC;
        }
302
303
304
305

        if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
            Vec3 boxVectors[3];
            force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
306
            gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
Rossen Apostolov's avatar
Rossen Apostolov committed
307
            gpuSetEwaldParameters(gpu);//, (float)alphaEwald, (int)kmax);
308
309
310
            method = EWALD;

        }
311
312
        data.nonbondedMethod = method;
        gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method);
313
314
315
316
    }

    // Initialize 1-4 nonbonded interactions.
    
317
    {
Peter Eastman's avatar
Peter Eastman committed
318
319
        vector<int> particle1(num14);
        vector<int> particle2(num14);
320
321
322
323
324
        vector<float> c6(num14);
        vector<float> c12(num14);
        vector<float> q1(num14);
        vector<float> q2(num14);
        for (int i = 0; i < num14; i++) {
325
            double charge, sig, eps;
326
            force.getExceptionParameters(nb14s[i], particle1[i], particle2[i], charge, sig, eps);
327
328
            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
329
330
            q1[i] = (float) charge;
            q2[i] = 1.0f;
331
        }
Peter Eastman's avatar
Peter Eastman committed
332
        gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, particle1, particle2, c6, c12, q1, q2);
333
334
335
    }
}

336
337
338
void CudaCalcNonbondedForceKernel::executeForces(OpenMMContextImpl& context) {
    if (data.primaryKernel == this)
        calcForces(context, data);
339
340
}

341
342
343
344
double CudaCalcNonbondedForceKernel::executeEnergy(OpenMMContextImpl& context) {
    if (data.primaryKernel == this)
        return calcEnergy(context, system);
    return 0.0;
345
346
}

347
CudaCalcGBSAOBCForceKernel::~CudaCalcGBSAOBCForceKernel() {
348
349
}

350
void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
351

Peter Eastman's avatar
Peter Eastman committed
352
    int numParticles = system.getNumParticles();
353
    _gpuContext* gpu = data.gpu;
Peter Eastman's avatar
Peter Eastman committed
354
355
    vector<float> radius(numParticles);
    vector<float> scale(numParticles);
356
    vector<float> charge(numParticles);
Peter Eastman's avatar
Peter Eastman committed
357
    for (int i = 0; i < numParticles; i++) {
358
359
        double particleCharge, particleRadius, scalingFactor;
        force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor);
Peter Eastman's avatar
Peter Eastman committed
360
        radius[i] = (float) particleRadius;
361
        scale[i] = (float) scalingFactor;
362
        charge[i] = (float) particleCharge;
363
    }
364
    gpuSetObcParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), radius, scale, charge);
365
366
}

367
void CudaCalcGBSAOBCForceKernel::executeForces(OpenMMContextImpl& context) {
368
369
}

370
static void initializeIntegration(const System& system, CudaPlatform::PlatformData& data, const Integrator& integrator) {
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387

    // 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>());
    }
388
389
390
    
    // Set masses.
    
Peter Eastman's avatar
Peter Eastman committed
391
392
393
394
    int numParticles = system.getNumParticles();
    vector<float> mass(numParticles);
    for (int i = 0; i < numParticles; i++)
        mass[i] = (float) system.getParticleMass(i);
395
396
397
398
    gpuSetMass(gpu, mass);
    
    // Set constraints.
    
399
    int numConstraints = system.getNumConstraints();
Peter Eastman's avatar
Peter Eastman committed
400
401
    vector<int> particle1(numConstraints);
    vector<int> particle2(numConstraints);
402
403
404
405
    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
406
        int particle1Index, particle2Index;
407
        double constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
408
409
410
        system.getConstraintParameters(i, particle1Index, particle2Index, constraintDistance);
        particle1[i] = particle1Index;
        particle2[i] = particle2Index;
411
        distance[i] = (float) constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
412
413
        invMass1[i] = 1.0f/mass[particle1Index];
        invMass2[i] = 1.0f/mass[particle2Index];
414
    }
415
    gpuSetConstraintParameters(gpu, particle1, particle2, distance, invMass1, invMass2, (float)integrator.getConstraintTolerance());
416
417
    
    // Finish initialization.
418

419
420
421
    gpuBuildThreadBlockWorkList(gpu);
    gpuBuildExclusionList(gpu);
    gpuBuildOutputBuffers(gpu);
422
    gpuSetConstants(gpu);
423
424
425
    kClearBornForces(gpu);
    kClearForces(gpu);
    cudaThreadSynchronize();
426
427
}

428
double CudaCalcGBSAOBCForceKernel::executeEnergy(OpenMMContextImpl& context) {
429
	return 0.0;
430
431
432
433
434
435
}

CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
}

void CudaIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
436
    initializeIntegration(system, data, integrator);
437
438
439
440
441
442
443
444
445
    prevStepSize = -1.0;
}

void CudaIntegrateVerletStepKernel::execute(OpenMMContextImpl& context, const VerletIntegrator& integrator) {
    _gpuContext* gpu = data.gpu;
    double stepSize = integrator.getStepSize();
    if (stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
446
        gpuSetVerletIntegrationParameters(gpu, (float) stepSize);
447
448
449
450
451
452
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevStepSize = stepSize;
    }
    kVerletUpdatePart1(gpu);
    kApplyFirstShake(gpu);
453
    kApplyFirstSettle(gpu);
454
    kApplyFirstCCMA(gpu);
455
    if (data.removeCM) {
456
        int step = (int) (context.getTime()/stepSize);
457
458
459
460
461
462
463
464
465
466
        if (step%data.cmMotionFrequency == 0)
            gpu->bCalculateCM = true;
    }
    kVerletUpdatePart2(gpu);
}

CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
}

void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
467
    initializeIntegration(system, data, integrator);
468
469
470
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
471
472
473
    prevStepSize = -1.0;
}

474
void CudaIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, const LangevinIntegrator& integrator) {
475
    _gpuContext* gpu = data.gpu;
476
477
478
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
479
480
481
482
    if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
        double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
483
        gpuSetIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature);
484
485
486
487
488
489
490
491
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
    kUpdatePart1(gpu);
    kApplyFirstShake(gpu);
492
    kApplyFirstSettle(gpu);
493
    kApplyFirstCCMA(gpu);
494
    if (data.removeCM) {
495
        int step = (int) (context.getTime()/stepSize);
496
497
498
        if (step%data.cmMotionFrequency == 0)
            gpu->bCalculateCM = true;
    }
499
500
    kUpdatePart2(gpu);
    kApplySecondShake(gpu);
501
    kApplySecondSettle(gpu);
502
    kApplySecondCCMA(gpu);
503
}
504
505
506
507
508

CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
}

void CudaIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
509
    initializeIntegration(system, data, integrator);
510
511
512
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
513
514
515
516
517
518
519
520
521
522
523
524
    prevStepSize = -1.0;
}

void CudaIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, const BrownianIntegrator& integrator) {
    _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);
525
        gpuSetBrownianIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature);
526
527
528
529
530
531
532
533
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
    kBrownianUpdatePart1(gpu);
    kApplyFirstShake(gpu);
534
    kApplyFirstSettle(gpu);
535
    kApplyFirstCCMA(gpu);
536
    if (data.removeCM) {
537
        int step = (int) (context.getTime()/stepSize);
538
539
540
541
542
        if (step%data.cmMotionFrequency == 0)
            gpu->bCalculateCM = true;
    }
    kBrownianUpdatePart2(gpu);
}
543
544
545
546
547

CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
}

void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
548
549
550
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) thermostat.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
551
552
553
554
555
    prevStepSize = -1.0;
}

void CudaApplyAndersenThermostatKernel::execute(OpenMMContextImpl& context) {
    _gpuContext* gpu = data.gpu;
556
557
    double temperature = context.getParameter(AndersenThermostat::Temperature());
    double frequency = context.getParameter(AndersenThermostat::CollisionFrequency());
558
559
560
561
    double stepSize = context.getIntegrator().getStepSize();
    if (temperature != prevTemp || frequency != prevFrequency || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
562
        gpuSetAndersenThermostatParameters(gpu, (float) temperature, (float) (frequency*stepSize));
563
564
565
566
567
568
569
570
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFrequency = frequency;
        prevStepSize = stepSize;
    }
    kCalculateAndersenThermostat(gpu);
}
571

572
void CudaCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
573
574
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
575
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
576
        masses[i] = system.getParticleMass(i);
577
578
}

579
double CudaCalcKineticEnergyKernel::execute(OpenMMContextImpl& context) {
580
581
582
    // We don't currently have a GPU kernel to do this, so we retrieve the velocities and calculate the energy
    // on the CPU.
    
583
    const Stream& velocities = context.getVelocities();
584
585
586
587
588
589
590
591
    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;
}
592

593
void CudaRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
594
    data.removeCM = true;
595
    data.cmMotionFrequency = force.getFrequency();
596
597
}

598
void CudaRemoveCMMotionKernel::execute(OpenMMContextImpl& context) {
599
}