CudaKernels.cpp 43 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
 * -------------------------------------------------------------------------- */

#include "CudaKernels.h"
28
#include "openmm/LangevinIntegrator.h"
29
#include "openmm/Context.h"
30
#include "openmm/OpenMMException.h"
31
#include "openmm/internal/ContextImpl.h"
32
#include "openmm/internal/NonbondedForceImpl.h"
33
#include "kernels/gputypes.h"
34
#include "kernels/cudaKernels.h"
35
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
36
37
#include <cmath>

38
extern "C" int OPENMMCUDA_EXPORT gpuSetConstants( gpuContext gpu );
39
40
41
42

using namespace OpenMM;
using namespace std;

43
44
45
46
void CudaCalcForcesAndEnergyKernel::initialize(const System& system) {
}

void CudaCalcForcesAndEnergyKernel::beginForceComputation(ContextImpl& context) {
47
    _gpuContext* gpu = data.gpu;
48
    if (data.nonbondedMethod != NO_CUTOFF && data.computeForceCount%100 == 0)
49
        gpuReorderAtoms(gpu);
50
    data.computeForceCount++;
51
52
53
54
    if (gpu->bIncludeGBSA || gpu->bIncludeGBVI)
        kClearBornSumAndForces(gpu);
    else
        kClearForces(gpu);
55
56
57
58
}

void CudaCalcForcesAndEnergyKernel::finishForceComputation(ContextImpl& context) {
    _gpuContext* gpu = data.gpu;
Mark Friedrichs's avatar
Mark Friedrichs committed
59
60

    if (gpu->bIncludeGBSA || gpu->bIncludeGBVI) {
61
        gpu->bRecalculateBornRadii = true;
62
63
        kCalculateCDLJObcGbsaForces1(gpu);
        kReduceObcGbsaBornForces(gpu);
Mark Friedrichs's avatar
Mark Friedrichs committed
64
65
66
67
68
        if (gpu->bIncludeGBSA ) { 
           kCalculateObcGbsaForces2(gpu);
        } else {
           kCalculateGBVIForces2(gpu);
        }
69
    }
70
    else if (data.hasNonbonded)
71
        kCalculateCDLJForces(gpu);
72
73
    if (data.hasCustomNonbonded)
        kCalculateCustomNonbondedForces(gpu, data.hasNonbonded);
74
    kCalculateLocalForces(gpu);
75
    kReduceForces(gpu);
76
77
}

78
void CudaCalcForcesAndEnergyKernel::beginEnergyComputation(ContextImpl& context) {
79
    _gpuContext* gpu = data.gpu;
80
81
82
83
    if (data.nonbondedMethod != NO_CUTOFF && data.stepCount%100 == 0)
        gpuReorderAtoms(gpu);
    data.stepCount++;
    kClearEnergy(gpu);
84
85
    if (gpu->bIncludeGBSA || gpu->bIncludeGBVI)
        kClearBornSumAndForces(gpu);
86
87
88
89
}

double CudaCalcForcesAndEnergyKernel::finishEnergyComputation(ContextImpl& context) {
    _gpuContext* gpu = data.gpu;
Mark Friedrichs's avatar
Mark Friedrichs committed
90
    if (gpu->bIncludeGBSA || gpu->bIncludeGBVI) {
91
92
93
        gpu->bRecalculateBornRadii = true;
        kCalculateCDLJObcGbsaForces1(gpu);
        kReduceObcGbsaBornForces(gpu);
Mark Friedrichs's avatar
Mark Friedrichs committed
94
95
96
97
98
        if (gpu->bIncludeGBSA ) {
           kCalculateObcGbsaForces2(gpu);
        } else {
           kCalculateGBVIForces2(gpu);
        }
99
    }
100
101
102
103
104
105
    else if (data.hasNonbonded)
        kCalculateCDLJForces(gpu);
    if (data.hasCustomNonbonded)
        kCalculateCustomNonbondedForces(gpu, data.hasNonbonded);
    kCalculateLocalForces(gpu);
    return kReduceEnergy(gpu)+data.ewaldSelfEnergy;
106
107
}

108
void CudaUpdateStateDataKernel::initialize(const System& system) {
109
110
}

111
double CudaUpdateStateDataKernel::getTime(const ContextImpl& context) const {
112
113
114
    return data.time;
}

115
void CudaUpdateStateDataKernel::setTime(ContextImpl& context, double time) {
116
117
118
    data.time = time;
}

119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
void CudaUpdateStateDataKernel::getPositions(ContextImpl& context, std::vector<Vec3>& positions) {
    _gpuContext* gpu = data.gpu;
    gpu->psPosq4->Download();
    int* order = gpu->psAtomIndex->_pSysData;
    int numParticles = context.getSystem().getNumParticles();
    positions.resize(numParticles);
    for (int i = 0; i < numParticles; ++i) {
        float4 pos = (*gpu->psPosq4)[i];
        int3 offset = gpu->posCellOffsets[i];
        positions[order[i]] = Vec3(pos.x-offset.x*gpu->sim.periodicBoxSizeX, pos.y-offset.y*gpu->sim.periodicBoxSizeY, pos.z-offset.z*gpu->sim.periodicBoxSizeZ);
    }
}

void CudaUpdateStateDataKernel::setPositions(ContextImpl& context, const std::vector<Vec3>& positions) {
    _gpuContext* gpu = data.gpu;
    int* order = gpu->psAtomIndex->_pSysData;
    int numParticles = context.getSystem().getNumParticles();
    for (int i = 0; i < numParticles; ++i) {
        float4& pos = (*gpu->psPosq4)[i];
        const Vec3& p = positions[order[i]];
139
140
141
        pos.x = (float) p[0];
        pos.y = (float) p[1];
        pos.z = (float) p[2];
142
143
    }
    gpu->psPosq4->Upload();
144
    for (int i = 0; i < (int) gpu->posCellOffsets.size(); i++)
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
        gpu->posCellOffsets[i] = make_int3(0, 0, 0);
}

void CudaUpdateStateDataKernel::getVelocities(ContextImpl& context, std::vector<Vec3>& velocities) {
    _gpuContext* gpu = data.gpu;
    gpu->psVelm4->Download();
    int* order = gpu->psAtomIndex->_pSysData;
    int numParticles = context.getSystem().getNumParticles();
    velocities.resize(numParticles);
    for (int i = 0; i < numParticles; ++i) {
        float4 vel = (*gpu->psVelm4)[i];
        velocities[order[i]] = Vec3(vel.x, vel.y, vel.z);
    }
}

void CudaUpdateStateDataKernel::setVelocities(ContextImpl& context, const std::vector<Vec3>& velocities) {
    _gpuContext* gpu = data.gpu;
    int* order = gpu->psAtomIndex->_pSysData;
    int numParticles = context.getSystem().getNumParticles();
    for (int i = 0; i < numParticles; ++i) {
        float4& vel = (*gpu->psVelm4)[i];
        const Vec3& v = velocities[order[i]];
167
168
169
        vel.x = (float) v[0];
        vel.y = (float) v[1];
        vel.z = (float) v[2];
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
    }
    gpu->psVelm4->Upload();
}

void CudaUpdateStateDataKernel::getForces(ContextImpl& context, std::vector<Vec3>& forces) {
    _gpuContext* gpu = data.gpu;
    int* order = gpu->psAtomIndex->_pSysData;
    gpu->psForce4->Download();
    int numParticles = context.getSystem().getNumParticles();
    forces.resize(numParticles);
    for (int i = 0; i < numParticles; ++i) {
        float4 force = (*gpu->psForce4)[i];
        forces[order[i]] = Vec3(force.x, force.y, force.z);
    }
}

186
187
188
189
190
191
192
193
194
195
196
197
198
void CudaUpdateStateDataKernel::getPeriodicBoxVectors(ContextImpl& context, Vec3& a, Vec3& b, Vec3& c) const {
    _gpuContext* gpu = data.gpu;
    a = Vec3(gpu->sim.periodicBoxSizeX, 0, 0);
    b = Vec3(0, gpu->sim.periodicBoxSizeY, 0);
    c = Vec3(0, 0, gpu->sim.periodicBoxSizeZ);
}

void CudaUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c) const {
    _gpuContext* gpu = data.gpu;
    gpuSetPeriodicBoxSize(gpu, a[0], b[1], c[2]);
    gpuSetConstants(gpu);
}

199
200
201
202
203
CudaCalcHarmonicBondForceKernel::~CudaCalcHarmonicBondForceKernel() {
}

void CudaCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
    data.hasBonds = true;
204
    numBonds = force.getNumBonds();
Peter Eastman's avatar
Peter Eastman committed
205
206
    vector<int> particle1(numBonds);
    vector<int> particle2(numBonds);
207
208
209
210
    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
211
        force.getBondParameters(i, particle1[i], particle2[i], lengthValue, kValue);
212
213
214
        length[i] = (float) lengthValue;
        k[i] = (float) kValue;
    }
Peter Eastman's avatar
Peter Eastman committed
215
    gpuSetBondParameters(data.gpu, particle1, particle2, length, k);
216
217
}

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

221
double CudaCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
222
223
224
    return 0.0;
}

225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
CudaCalcCustomBondForceKernel::~CudaCalcCustomBondForceKernel() {
}

void CudaCalcCustomBondForceKernel::initialize(const System& system, const CustomBondForce& force) {
    numBonds = force.getNumBonds();
    vector<int> particle1(numBonds);
    vector<int> particle2(numBonds);
    vector<vector<double> > params(numBonds);
    for (int i = 0; i < numBonds; i++)
        force.getBondParameters(i, particle1[i], particle2[i], params[i]);
    vector<string> paramNames;
    for (int i = 0; i < force.getNumPerBondParameters(); i++)
        paramNames.push_back(force.getPerBondParameterName(i));
    globalParamNames.resize(force.getNumGlobalParameters());
    globalParamValues.resize(force.getNumGlobalParameters());
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
        globalParamNames[i] = force.getGlobalParameterName(i);
        globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
    }
    gpuSetCustomBondParameters(data.gpu, particle1, particle2, params, force.getEnergyFunction(), paramNames, globalParamNames);
    if (globalParamValues.size() > 0)
Peter Eastman's avatar
Peter Eastman committed
246
        SetCustomBondGlobalParams(globalParamValues);
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
}

void CudaCalcCustomBondForceKernel::executeForces(ContextImpl& context) {
    updateGlobalParams(context);
    kCalculateCustomBondForces(data.gpu);
}

double CudaCalcCustomBondForceKernel::executeEnergy(ContextImpl& context) {
    updateGlobalParams(context);
    kCalculateCustomBondForces(data.gpu);
    return 0.0;
}

void CudaCalcCustomBondForceKernel::updateGlobalParams(ContextImpl& context) {
    bool changed = false;
    for (int i = 0; i < (int) globalParamNames.size(); i++) {
        float value = (float) context.getParameter(globalParamNames[i]);
        if (value != globalParamValues[i])
            changed = true;
        globalParamValues[i] = value;
    }
    if (changed)
Peter Eastman's avatar
Peter Eastman committed
269
        SetCustomBondGlobalParams(globalParamValues);
270
271
}

272
273
274
275
276
CudaCalcHarmonicAngleForceKernel::~CudaCalcHarmonicAngleForceKernel() {
}

void CudaCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
    data.hasAngles = true;
277
    numAngles = force.getNumAngles();
278
    const float RadiansToDegrees = (float) (180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
279
280
281
    vector<int> particle1(numAngles);
    vector<int> particle2(numAngles);
    vector<int> particle3(numAngles);
282
283
284
285
    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
286
        force.getAngleParameters(i, particle1[i], particle2[i], particle3[i], angleValue, kValue);
287
288
        angle[i] = (float) (angleValue*RadiansToDegrees);
        k[i] = (float) kValue;
289
    }
Peter Eastman's avatar
Peter Eastman committed
290
    gpuSetBondAngleParameters(data.gpu, particle1, particle2, particle3, angle, k);
291
}
292

293
void CudaCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
294
295
}

296
double CudaCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
297
298
299
    return 0.0;
}

300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
CudaCalcCustomAngleForceKernel::~CudaCalcCustomAngleForceKernel() {
}

void CudaCalcCustomAngleForceKernel::initialize(const System& system, const CustomAngleForce& force) {
    numAngles = force.getNumAngles();
    vector<int> particle1(numAngles);
    vector<int> particle2(numAngles);
    vector<int> particle3(numAngles);
    vector<vector<double> > params(numAngles);
    for (int i = 0; i < numAngles; i++)
        force.getAngleParameters(i, particle1[i], particle2[i], particle3[i], params[i]);
    vector<string> paramNames;
    for (int i = 0; i < force.getNumPerAngleParameters(); i++)
        paramNames.push_back(force.getPerAngleParameterName(i));
    globalParamNames.resize(force.getNumGlobalParameters());
    globalParamValues.resize(force.getNumGlobalParameters());
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
        globalParamNames[i] = force.getGlobalParameterName(i);
        globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
    }
    gpuSetCustomAngleParameters(data.gpu, particle1, particle2, particle3, params, force.getEnergyFunction(), paramNames, globalParamNames);
    if (globalParamValues.size() > 0)
        SetCustomAngleGlobalParams(globalParamValues);
}

void CudaCalcCustomAngleForceKernel::executeForces(ContextImpl& context) {
    updateGlobalParams(context);
    kCalculateCustomAngleForces(data.gpu);
}

double CudaCalcCustomAngleForceKernel::executeEnergy(ContextImpl& context) {
    updateGlobalParams(context);
    kCalculateCustomAngleForces(data.gpu);
    return 0.0;
}

void CudaCalcCustomAngleForceKernel::updateGlobalParams(ContextImpl& context) {
    bool changed = false;
    for (int i = 0; i < (int) globalParamNames.size(); i++) {
        float value = (float) context.getParameter(globalParamNames[i]);
        if (value != globalParamValues[i])
            changed = true;
        globalParamValues[i] = value;
    }
    if (changed)
        SetCustomAngleGlobalParams(globalParamValues);
}

348
349
350
351
352
353
CudaCalcPeriodicTorsionForceKernel::~CudaCalcPeriodicTorsionForceKernel() {
}

void CudaCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
    data.hasPeriodicTorsions = true;
    numTorsions = force.getNumTorsions();
354
    const float RadiansToDegrees = (float)(180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
355
356
357
358
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
359
360
361
362
363
    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
364
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], periodicity[i], phaseValue, kValue);
365
366
        k[i] = (float) kValue;
        phase[i] = (float) (phaseValue*RadiansToDegrees);
367
    }
Peter Eastman's avatar
Peter Eastman committed
368
    gpuSetDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, k, phase, periodicity);
369
370
}

371
void CudaCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
372
373
}

374
double CudaCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
375
376
377
378
379
380
381
382
383
    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
384
385
386
387
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
388
389
390
391
392
393
394
395
    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
396
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], c[0], c[1], c[2], c[3], c[4], c[5]);
397
398
399
400
401
402
        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];
403
    }
Peter Eastman's avatar
Peter Eastman committed
404
    gpuSetRbDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
405
406
}

407
void CudaCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
408
409
}

410
double CudaCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
411
412
413
    return 0.0;
}

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
CudaCalcCustomTorsionForceKernel::~CudaCalcCustomTorsionForceKernel() {
}

void CudaCalcCustomTorsionForceKernel::initialize(const System& system, const CustomTorsionForce& force) {
    numTorsions = force.getNumTorsions();
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
    vector<vector<double> > params(numTorsions);
    for (int i = 0; i < numTorsions; i++)
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], params[i]);
    vector<string> paramNames;
    for (int i = 0; i < force.getNumPerTorsionParameters(); i++)
        paramNames.push_back(force.getPerTorsionParameterName(i));
    globalParamNames.resize(force.getNumGlobalParameters());
    globalParamValues.resize(force.getNumGlobalParameters());
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
        globalParamNames[i] = force.getGlobalParameterName(i);
        globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
    }
    gpuSetCustomTorsionParameters(data.gpu, particle1, particle2, particle3, particle4, params, force.getEnergyFunction(), paramNames, globalParamNames);
    if (globalParamValues.size() > 0)
        SetCustomTorsionGlobalParams(globalParamValues);
}

void CudaCalcCustomTorsionForceKernel::executeForces(ContextImpl& context) {
    updateGlobalParams(context);
    kCalculateCustomTorsionForces(data.gpu);
}

double CudaCalcCustomTorsionForceKernel::executeEnergy(ContextImpl& context) {
    updateGlobalParams(context);
    kCalculateCustomTorsionForces(data.gpu);
    return 0.0;
}

void CudaCalcCustomTorsionForceKernel::updateGlobalParams(ContextImpl& context) {
    bool changed = false;
    for (int i = 0; i < (int) globalParamNames.size(); i++) {
        float value = (float) context.getParameter(globalParamNames[i]);
        if (value != globalParamValues[i])
            changed = true;
        globalParamValues[i] = value;
    }
    if (changed)
        SetCustomTorsionGlobalParams(globalParamValues);
}

463
464
465
CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() {
}

466
void CudaCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
467
    data.hasNonbonded = true;
Peter Eastman's avatar
Peter Eastman committed
468
    numParticles = force.getNumParticles();
469
    _gpuContext* gpu = data.gpu;
470
471
472
473

    // Identify which exceptions are 1-4 interactions.

    vector<pair<int, int> > exclusions;
474
    vector<int> exceptions;
475
476
477
478
479
480
    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)
481
            exceptions.push_back(i);
482
483
    }

484
485
    // Initialize nonbonded interactions.
    
486
    {
Peter Eastman's avatar
Peter Eastman committed
487
488
489
490
        vector<int> particle(numParticles);
        vector<float> c6(numParticles);
        vector<float> c12(numParticles);
        vector<float> q(numParticles);
491
        vector<char> symbol;
Peter Eastman's avatar
Peter Eastman committed
492
493
        vector<vector<int> > exclusionList(numParticles);
        for (int i = 0; i < numParticles; i++) {
494
            double charge, radius, depth;
Peter Eastman's avatar
Peter Eastman committed
495
496
            force.getParticleParameters(i, charge, radius, depth);
            particle[i] = i;
497
498
499
            q[i] = (float) charge;
            c6[i] = (float) (4*depth*pow(radius, 6.0));
            c12[i] = (float) (4*depth*pow(radius, 12.0));
500
501
            exclusionList[i].push_back(i);
        }
502
        for (int i = 0; i < (int)exclusions.size(); i++) {
503
504
505
            exclusionList[exclusions[i].first].push_back(exclusions[i].second);
            exclusionList[exclusions[i].second].push_back(exclusions[i].first);
        }
506
507
        CudaNonbondedMethod method = NO_CUTOFF;
        if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
508
            gpuSetNonbondedCutoff(gpu, (float) force.getCutoffDistance(), (float) force.getReactionFieldDielectric());
509
510
511
512
513
            method = CUTOFF;
        }
        if (force.getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
            method = PERIODIC;
        }
514
515
        if (force.getNonbondedMethod() == NonbondedForce::Ewald || force.getNonbondedMethod() == NonbondedForce::PME) {
            if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
516
517
518
                double alpha;
                int kmaxx, kmaxy, kmaxz;
                NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmaxx, kmaxy, kmaxz);
519
520
521
522
                gpuSetEwaldParameters(gpu, (float) alpha, kmaxx, kmaxy, kmaxz);
                method = EWALD;
            }
            else {
523
524
525
                double alpha;
                int gridSizeX, gridSizeY, gridSizeZ;
                NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSizeX, gridSizeY, gridSizeZ);
Peter Eastman's avatar
Peter Eastman committed
526
                gpuSetPMEParameters(gpu, (float) alpha, gridSizeX, gridSizeY, gridSizeZ);
527
528
                method = PARTICLE_MESH_EWALD;
            }
529
        }
530
        data.nonbondedMethod = method;
531
        gpuSetCoulombParameters(gpu, (float) ONE_4PI_EPS0, particle, c6, c12, q, symbol, exclusionList, method);
532
533
534
535

        // Compute the Ewald self energy.

        data.ewaldSelfEnergy = 0.0;
536
537
538
539
540
        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];
        }
541
542
543
544
    }

    // Initialize 1-4 nonbonded interactions.
    
545
    {
546
547
548
549
550
551
552
553
        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++) {
554
            double charge, sig, eps;
555
            force.getExceptionParameters(exceptions[i], particle1[i], particle2[i], charge, sig, eps);
556
557
            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
558
559
            q1[i] = (float) charge;
            q2[i] = 1.0f;
560
        }
561
        gpuSetLJ14Parameters(gpu, (float) ONE_4PI_EPS0, 1.0f, particle1, particle2, c6, c12, q1, q2);
562
563
564
    }
}

565
void CudaCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
566
567
}

568
double CudaCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
569
    return 0.0;
570
571
}

572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
CudaCalcCustomNonbondedForceKernel::~CudaCalcCustomNonbondedForceKernel() {
}

void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
    data.hasCustomNonbonded = true;
    numParticles = force.getNumParticles();
    _gpuContext* gpu = data.gpu;

    // 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);
    }
590
591
592
593
594
    for (int i = 0; i < force.getNumExclusions(); i++) {
        int particle1, particle2;
        force.getExclusionParticles(i, particle1, particle2);
        exclusionList[particle1].push_back(particle2);
        exclusionList[particle2].push_back(particle1);
595
596
597
598
599
600
601
602
603
    }
    CudaNonbondedMethod method = NO_CUTOFF;
    if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
        method = CUTOFF;
    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic) {
        method = PERIODIC;
    }
    data.customNonbondedMethod = method;

604
605
606
607
608
609
610
611
612
613
614
    // 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);
    }

615
616
617
    // Record information for the expressions.

    vector<string> paramNames;
618
619
    for (int i = 0; i < force.getNumPerParticleParameters(); i++)
        paramNames.push_back(force.getPerParticleParameterName(i));
620
621
622
623
    globalParamNames.resize(force.getNumGlobalParameters());
    globalParamValues.resize(force.getNumGlobalParameters());
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
        globalParamNames[i] = force.getGlobalParameterName(i);
624
        globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
625
    }
626
    gpuSetCustomNonbondedParameters(gpu, parameters, exclusionList, method, (float) force.getCutoffDistance(), force.getEnergyFunction(), paramNames, globalParamNames);
627
    if (globalParamValues.size() > 0)
Peter Eastman's avatar
Peter Eastman committed
628
        SetCustomNonbondedGlobalParams(globalParamValues);
629
630
631
}

void CudaCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
632
    updateGlobalParams(context);
633
634
635
}

double CudaCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
636
    updateGlobalParams(context);
637
638
639
    return 0.0;
}

640
641
void CudaCalcCustomNonbondedForceKernel::updateGlobalParams(ContextImpl& context) {
    bool changed = false;
642
    for (int i = 0; i < (int) globalParamNames.size(); i++) {
643
644
645
646
647
648
        float value = (float) context.getParameter(globalParamNames[i]);
        if (value != globalParamValues[i])
            changed = true;
        globalParamValues[i] = value;
    }
    if (changed)
Peter Eastman's avatar
Peter Eastman committed
649
        SetCustomNonbondedGlobalParams(globalParamValues);
650
651
}

652
CudaCalcGBSAOBCForceKernel::~CudaCalcGBSAOBCForceKernel() {
653
654
}

655
void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
656

Peter Eastman's avatar
Peter Eastman committed
657
    int numParticles = system.getNumParticles();
658
    _gpuContext* gpu = data.gpu;
Peter Eastman's avatar
Peter Eastman committed
659
660
    vector<float> radius(numParticles);
    vector<float> scale(numParticles);
661
    vector<float> charge(numParticles);
Peter Eastman's avatar
Peter Eastman committed
662
    for (int i = 0; i < numParticles; i++) {
663
664
        double particleCharge, particleRadius, scalingFactor;
        force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor);
Peter Eastman's avatar
Peter Eastman committed
665
        radius[i] = (float) particleRadius;
666
        scale[i] = (float) scalingFactor;
667
        charge[i] = (float) particleCharge;
668
    }
669
    gpuSetObcParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), radius, scale, charge);
670
671
}

672
void CudaCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
673
674
}

Mark Friedrichs's avatar
Mark Friedrichs committed
675
676
677
678
679
680
681
682
683
684
685
686
687
688
CudaCalcGBVIForceKernel::~CudaCalcGBVIForceKernel() {
}

void CudaCalcGBVIForceKernel::initialize(const System& system, const GBVIForce& force, const std::vector<double> & inputScaledRadii) {

    int numParticles = system.getNumParticles();
    _gpuContext* gpu = data.gpu;

    vector<int> particle(numParticles);
    vector<float> radius(numParticles);
    vector<float> scaledRadii(numParticles);
    vector<float> gammas(numParticles);

    for (int i = 0; i < numParticles; i++) {
689
        double charge, particleRadius, gamma;
Mark Friedrichs's avatar
Mark Friedrichs committed
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
        force.getParticleParameters(i, charge, particleRadius, gamma );
        particle[i]                  = i;
        radius[i]                    = (float) particleRadius;
        gammas[i]                    = (float) gamma;
        scaledRadii[i]               = (float) inputScaledRadii[i];
    }
    gpuSetGBVIParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), particle,
                         radius, gammas, scaledRadii );
}

void CudaCalcGBVIForceKernel::executeForces(ContextImpl& context) {
}

double CudaCalcGBVIForceKernel::executeEnergy(ContextImpl& context) {
    return 0.0;
}

707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
CudaCalcCustomExternalForceKernel::~CudaCalcCustomExternalForceKernel() {
}

void CudaCalcCustomExternalForceKernel::initialize(const System& system, const CustomExternalForce& force) {
    numParticles = force.getNumParticles();
    vector<int> particle(numParticles);
    vector<vector<double> > params(numParticles);
    for (int i = 0; i < numParticles; i++)
        force.getParticleParameters(i, particle[i], params[i]);
    vector<string> paramNames;
    for (int i = 0; i < force.getNumPerParticleParameters(); i++)
        paramNames.push_back(force.getPerParticleParameterName(i));
    globalParamNames.resize(force.getNumGlobalParameters());
    globalParamValues.resize(force.getNumGlobalParameters());
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
        globalParamNames[i] = force.getGlobalParameterName(i);
        globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
    }
    gpuSetCustomExternalParameters(data.gpu, particle, params, force.getEnergyFunction(), paramNames, globalParamNames);
    if (globalParamValues.size() > 0)
Peter Eastman's avatar
Peter Eastman committed
727
        SetCustomExternalGlobalParams(globalParamValues);
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
}

void CudaCalcCustomExternalForceKernel::executeForces(ContextImpl& context) {
    updateGlobalParams(context);
    kCalculateCustomExternalForces(data.gpu);
}

double CudaCalcCustomExternalForceKernel::executeEnergy(ContextImpl& context) {
    updateGlobalParams(context);
    kCalculateCustomExternalForces(data.gpu);
    return 0.0;
}

void CudaCalcCustomExternalForceKernel::updateGlobalParams(ContextImpl& context) {
    bool changed = false;
    for (int i = 0; i < (int) globalParamNames.size(); i++) {
        float value = (float) context.getParameter(globalParamNames[i]);
        if (value != globalParamValues[i])
            changed = true;
        globalParamValues[i] = value;
    }
    if (changed)
Peter Eastman's avatar
Peter Eastman committed
750
        SetCustomExternalGlobalParams(globalParamValues);
751
752
}

753
void OPENMMCUDA_EXPORT OpenMM::cudaOpenMMInitializeIntegration(const System& system, CudaPlatform::PlatformData& data, const Integrator& integrator) {
754
755
756
757
758
759
760
761
762
763
764
765
766
767

    // 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) {
768
769
        gpuSetCoulombParameters(gpu, (float) ONE_4PI_EPS0, vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<char>(), vector<vector<int> >(), NO_CUTOFF);
        gpuSetLJ14Parameters(gpu, (float) ONE_4PI_EPS0, 1.0f, vector<int>(), vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<float>());
Mark Friedrichs's avatar
Update  
Mark Friedrichs committed
770
/*
771
772
        if (gpu->bIncludeGBSA || gpu->bIncludeGBVI)
            throw OpenMMException("CudaPlatform requires GBSAOBCForce and GBVIForce to be used with a NonbondedForce");
Mark Friedrichs's avatar
Update  
Mark Friedrichs committed
773
*/
774
    }
775
776
777
    
    // Set masses.
    
Peter Eastman's avatar
Peter Eastman committed
778
779
780
781
    int numParticles = system.getNumParticles();
    vector<float> mass(numParticles);
    for (int i = 0; i < numParticles; i++)
        mass[i] = (float) system.getParticleMass(i);
782
783
784
785
    gpuSetMass(gpu, mass);
    
    // Set constraints.
    
786
    int numConstraints = system.getNumConstraints();
Peter Eastman's avatar
Peter Eastman committed
787
788
    vector<int> particle1(numConstraints);
    vector<int> particle2(numConstraints);
789
790
791
792
    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
793
        int particle1Index, particle2Index;
794
        double constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
795
796
797
        system.getConstraintParameters(i, particle1Index, particle2Index, constraintDistance);
        particle1[i] = particle1Index;
        particle2[i] = particle2Index;
798
        distance[i] = (float) constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
799
800
        invMass1[i] = 1.0f/mass[particle1Index];
        invMass2[i] = 1.0f/mass[particle2Index];
801
    }
802
    gpuSetConstraintParameters(gpu, particle1, particle2, distance, invMass1, invMass2, (float)integrator.getConstraintTolerance());
803
804
    
    // Finish initialization.
805

806
807
808
    gpuBuildThreadBlockWorkList(gpu);
    gpuBuildExclusionList(gpu);
    gpuBuildOutputBuffers(gpu);
809
    gpuSetConstants(gpu);
810
811
812
813
    if (gpu->bIncludeGBSA || gpu->bIncludeGBVI)
        kClearBornSumAndForces(gpu);
    else
        kClearForces(gpu);
814
    cudaThreadSynchronize();
815
816
}

817
double CudaCalcGBSAOBCForceKernel::executeEnergy(ContextImpl& context) {
818
	return 0.0;
819
820
821
822
823
824
}

CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
}

void CudaIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
825
    cudaOpenMMInitializeIntegration(system, data, integrator);
826
827
828
    prevStepSize = -1.0;
}

829
void CudaIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
830
831
832
833
834
    _gpuContext* gpu = data.gpu;
    double stepSize = integrator.getStepSize();
    if (stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
835
        gpuSetVerletIntegrationParameters(gpu, (float) stepSize, 0.0f);
836
837
838
839
840
        gpuSetConstants(gpu);
        prevStepSize = stepSize;
    }
    kVerletUpdatePart1(gpu);
    kApplyFirstShake(gpu);
841
    kApplyFirstSettle(gpu);
842
    kApplyFirstCCMA(gpu);
843
844
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
845
846
            gpu->bCalculateCM = true;
    kVerletUpdatePart2(gpu);
847
    data.time += stepSize;
848
    data.stepCount++;
849
850
851
852
853
854
}

CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
}

void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
855
    cudaOpenMMInitializeIntegration(system, data, integrator);
856
857
858
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
859
860
    prevTemp = -1.0;
    prevFriction = -1.0;
861
862
863
    prevStepSize = -1.0;
}

864
void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
865
    _gpuContext* gpu = data.gpu;
866
867
868
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
869
870
    if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
871
        
872
        double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
873
        gpuSetLangevinIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature, 0.0f);
874
875
876
877
878
879
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
880
    kLangevinUpdatePart1(gpu);
881
882
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
883
            gpu->bCalculateCM = true;
884
    kLangevinUpdatePart2(gpu);
885
    kApplySecondShake(gpu);
886
    kApplySecondSettle(gpu);
887
    kApplySecondCCMA(gpu);
888
    kSetVelocitiesFromPositions(gpu);
889
    data.time += stepSize;
890
    data.stepCount++;
891
}
892
893
894
895
896

CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
}

void CudaIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
897
    cudaOpenMMInitializeIntegration(system, data, integrator);
898
899
900
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
901
902
    prevTemp = -1.0;
    prevFriction = -1.0;
903
904
905
    prevStepSize = -1.0;
}

906
void CudaIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
907
908
909
910
911
912
913
914
    _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);
915
        gpuSetBrownianIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature);
916
917
918
919
920
921
922
923
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
    kBrownianUpdatePart1(gpu);
    kApplyFirstShake(gpu);
924
    kApplyFirstSettle(gpu);
925
    kApplyFirstCCMA(gpu);
926
927
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
928
929
            gpu->bCalculateCM = true;
    kBrownianUpdatePart2(gpu);
930
    data.time += stepSize;
931
932
933
934
935
936
937
    data.stepCount++;
}

CudaIntegrateVariableVerletStepKernel::~CudaIntegrateVariableVerletStepKernel() {
}

void CudaIntegrateVariableVerletStepKernel::initialize(const System& system, const VariableVerletIntegrator& integrator) {
938
    cudaOpenMMInitializeIntegration(system, data, integrator);
939
940
941
    prevErrorTol = -1.0;
}

942
void CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
943
944
945
946
947
948
949
950
951
    _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;
    }
952
    float maxStepSize = (float)(maxTime-data.time);
953
    kSelectVerletStepSize(gpu, maxStepSize);
954
955
956
957
958
959
960
961
962
963
    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;
964
965
    if ((*gpu->psStepSize)[0].y == maxStepSize)
        data.time = maxTime; // Avoid round-off error
966
    data.stepCount++;
967
}
968

969
970
971
972
CudaIntegrateVariableLangevinStepKernel::~CudaIntegrateVariableLangevinStepKernel() {
}

void CudaIntegrateVariableLangevinStepKernel::initialize(const System& system, const VariableLangevinIntegrator& integrator) {
973
    cudaOpenMMInitializeIntegration(system, data, integrator);
974
975
976
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
977
978
    prevTemp = -1.0;
    prevFriction = -1.0;
979
980
981
982
983
984
985
986
987
988
989
990
    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);
991
        gpuSetLangevinIntegrationParameters(gpu, (float) tau, 0.0f, (float) temperature, (float) errorTol);
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevErrorTol = errorTol;
    }
    float maxStepSize = (float)(maxTime-data.time);
    kSelectLangevinStepSize(gpu, maxStepSize);
    kLangevinUpdatePart1(gpu);
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
            gpu->bCalculateCM = true;
    kLangevinUpdatePart2(gpu);
    kApplySecondShake(gpu);
    kApplySecondSettle(gpu);
    kApplySecondCCMA(gpu);
1008
    kSetVelocitiesFromPositions(gpu);
1009
1010
1011
1012
1013
1014
1015
    gpu->psStepSize->Download();
    data.time += (*gpu->psStepSize)[0].y;
    if ((*gpu->psStepSize)[0].y == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
}

1016
1017
1018
1019
CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
}

void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
1020
1021
1022
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) thermostat.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
1023
1024
    prevTemp = -1.0;
    prevFrequency = -1.0;
1025
1026
1027
    prevStepSize = -1.0;
}

1028
void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
1029
    _gpuContext* gpu = data.gpu;
1030
1031
    double temperature = context.getParameter(AndersenThermostat::Temperature());
    double frequency = context.getParameter(AndersenThermostat::CollisionFrequency());
1032
1033
1034
1035
    double stepSize = context.getIntegrator().getStepSize();
    if (temperature != prevTemp || frequency != prevFrequency || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
1036
        gpuSetAndersenThermostatParameters(gpu, (float) temperature, (float) frequency);
1037
1038
1039
1040
1041
1042
1043
1044
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFrequency = frequency;
        prevStepSize = stepSize;
    }
    kCalculateAndersenThermostat(gpu);
}
1045

1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
CudaApplyMonteCarloBarostatKernel::~CudaApplyMonteCarloBarostatKernel() {
    if (moleculeAtoms != NULL)
        delete moleculeAtoms;
    if (moleculeStartIndex != NULL)
        delete moleculeStartIndex;
}

void CudaApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& thermostat) {
}

void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
    if (!hasInitializedMolecules) {
        hasInitializedMolecules = true;

        // Create the arrays with the molecule definitions.

        vector<vector<int> > molecules = context.getMolecules();
        numMolecules = molecules.size();
        moleculeAtoms = new CUDAStream<int>(context.getSystem().getNumParticles(), 1, "moleculeAtoms");
        moleculeStartIndex = new CUDAStream<int>(numMolecules+1, 1, "moleculeStartIndex");
        int index = 0;
        for (int i = 0; i < numMolecules; i++) {
            (*moleculeStartIndex)[i] = index;
            for (int j = 0; j < (int) molecules[i].size(); j++)
                (*moleculeAtoms)[index++] = molecules[i][j];
        }
        (*moleculeStartIndex)[numMolecules] = index;
        moleculeAtoms->Upload();
        moleculeStartIndex->Upload();
    }
    _gpuContext* gpu = data.gpu;
    gpu->psPosqP4->CopyFrom(*gpu->psPosq4);
    kScaleAtomCoordinates(gpu, scale, *moleculeAtoms, *moleculeStartIndex);
}

void CudaApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
    _gpuContext* gpu = data.gpu;
    gpu->psPosq4->CopyFrom(*gpu->psPosqP4);
}

1086
void CudaCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
1087
1088
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
1089
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1090
        masses[i] = system.getParticleMass(i);
1091
1092
}

1093
double CudaCalcKineticEnergyKernel::execute(ContextImpl& context) {
1094
1095
1096
    // We don't currently have a GPU kernel to do this, so we retrieve the velocities and calculate the energy
    // on the CPU.
    
1097
1098
    _gpuContext* gpu = data.gpu;
    gpu->psVelm4->Download();
1099
    double energy = 0.0;
1100
1101
1102
1103
    for (int i = 0; i < (int) masses.size(); ++i) {
        float4 v = (*gpu->psVelm4)[i];
        energy += masses[i]*(v.x*v.x+v.y*v.y+v.z*v.z);
    }
1104
1105
    return 0.5*energy;
}
1106

1107
void CudaRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
1108
    data.removeCM = true;
1109
    data.cmMotionFrequency = force.getFrequency();
1110
1111
}

1112
void CudaRemoveCMMotionKernel::execute(ContextImpl& context) {
1113
}