CudaKernels.cpp 43.8 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
    else if (data.hasNonbonded)
        kCalculateCDLJForces(gpu);
    if (data.hasCustomNonbonded)
        kCalculateCustomNonbondedForces(gpu, data.hasNonbonded);
    kCalculateLocalForces(gpu);
105
106
107
108
    double energy = kReduceEnergy(gpu)+data.ewaldSelfEnergy;
    if (data.dispersionCoefficient != 0.0)
        energy += data.dispersionCoefficient/(gpu->sim.periodicBoxSizeX*gpu->sim.periodicBoxSizeY*gpu->sim.periodicBoxSizeZ);
    return energy;
109
110
}

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

114
double CudaUpdateStateDataKernel::getTime(const ContextImpl& context) const {
115
116
117
    return data.time;
}

118
void CudaUpdateStateDataKernel::setTime(ContextImpl& context, double time) {
119
120
121
    data.time = time;
}

122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
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]];
142
143
144
        pos.x = (float) p[0];
        pos.y = (float) p[1];
        pos.z = (float) p[2];
145
146
    }
    gpu->psPosq4->Upload();
147
    for (int i = 0; i < (int) gpu->posCellOffsets.size(); i++)
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
        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]];
170
171
172
        vel.x = (float) v[0];
        vel.y = (float) v[1];
        vel.z = (float) v[2];
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
    }
    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);
    }
}

189
190
191
192
193
194
195
196
197
198
199
200
201
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);
}

202
203
204
205
206
207
208
void CudaApplyConstraintsKernel::initialize(const System& system) {
}

void CudaApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
    kApplyConstraints(data.gpu);
}

209
210
211
212
213
CudaCalcHarmonicBondForceKernel::~CudaCalcHarmonicBondForceKernel() {
}

void CudaCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
    data.hasBonds = true;
214
    numBonds = force.getNumBonds();
Peter Eastman's avatar
Peter Eastman committed
215
216
    vector<int> particle1(numBonds);
    vector<int> particle2(numBonds);
217
218
219
220
    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
221
        force.getBondParameters(i, particle1[i], particle2[i], lengthValue, kValue);
222
223
224
        length[i] = (float) lengthValue;
        k[i] = (float) kValue;
    }
Peter Eastman's avatar
Peter Eastman committed
225
    gpuSetBondParameters(data.gpu, particle1, particle2, length, k);
226
227
}

228
void CudaCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
229
230
}

231
double CudaCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
232
233
234
    return 0.0;
}

235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
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
256
        SetCustomBondGlobalParams(globalParamValues);
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
}

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
279
        SetCustomBondGlobalParams(globalParamValues);
280
281
}

282
283
284
285
286
CudaCalcHarmonicAngleForceKernel::~CudaCalcHarmonicAngleForceKernel() {
}

void CudaCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
    data.hasAngles = true;
287
    numAngles = force.getNumAngles();
288
    const float RadiansToDegrees = (float) (180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
289
290
291
    vector<int> particle1(numAngles);
    vector<int> particle2(numAngles);
    vector<int> particle3(numAngles);
292
293
294
295
    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
296
        force.getAngleParameters(i, particle1[i], particle2[i], particle3[i], angleValue, kValue);
297
298
        angle[i] = (float) (angleValue*RadiansToDegrees);
        k[i] = (float) kValue;
299
    }
Peter Eastman's avatar
Peter Eastman committed
300
    gpuSetBondAngleParameters(data.gpu, particle1, particle2, particle3, angle, k);
301
}
302

303
void CudaCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
304
305
}

306
double CudaCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
307
308
309
    return 0.0;
}

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
348
349
350
351
352
353
354
355
356
357
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);
}

358
359
360
361
362
363
CudaCalcPeriodicTorsionForceKernel::~CudaCalcPeriodicTorsionForceKernel() {
}

void CudaCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
    data.hasPeriodicTorsions = true;
    numTorsions = force.getNumTorsions();
364
    const float RadiansToDegrees = (float)(180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
365
366
367
368
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
369
370
371
372
373
    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
374
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], periodicity[i], phaseValue, kValue);
375
376
        k[i] = (float) kValue;
        phase[i] = (float) (phaseValue*RadiansToDegrees);
377
    }
Peter Eastman's avatar
Peter Eastman committed
378
    gpuSetDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, k, phase, periodicity);
379
380
}

381
void CudaCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
382
383
}

384
double CudaCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
385
386
387
388
389
390
391
392
393
    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
394
395
396
397
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
398
399
400
401
402
403
404
405
    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
406
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], c[0], c[1], c[2], c[3], c[4], c[5]);
407
408
409
410
411
412
        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];
413
    }
Peter Eastman's avatar
Peter Eastman committed
414
    gpuSetRbDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
415
416
}

417
void CudaCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
418
419
}

420
double CudaCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
421
422
423
    return 0.0;
}

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
466
467
468
469
470
471
472
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);
}

473
474
475
CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() {
}

476
void CudaCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
477
    data.hasNonbonded = true;
Peter Eastman's avatar
Peter Eastman committed
478
    numParticles = force.getNumParticles();
479
    _gpuContext* gpu = data.gpu;
480
481
482
483

    // Identify which exceptions are 1-4 interactions.

    vector<pair<int, int> > exclusions;
484
    vector<int> exceptions;
485
486
487
488
489
490
    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)
491
            exceptions.push_back(i);
492
493
    }

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

        // Compute the Ewald self energy.

        data.ewaldSelfEnergy = 0.0;
546
547
548
549
550
        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];
        }
551
552
553
554
555
556
557

        // Compute the long range dispersion correction.

        if (force.getUseDispersionCorrection())
            data.dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
        else
            data.dispersionCoefficient = 0.0;
558
559
560
561
    }

    // Initialize 1-4 nonbonded interactions.
    
562
    {
563
564
565
566
567
568
569
570
        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++) {
571
            double charge, sig, eps;
572
            force.getExceptionParameters(exceptions[i], particle1[i], particle2[i], charge, sig, eps);
573
574
            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
575
576
            q1[i] = (float) charge;
            q2[i] = 1.0f;
577
        }
578
        gpuSetLJ14Parameters(gpu, (float) ONE_4PI_EPS0, 1.0f, particle1, particle2, c6, c12, q1, q2);
579
580
581
    }
}

582
void CudaCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
583
584
}

585
double CudaCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
586
    return 0.0;
587
588
}

589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
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);
    }
607
608
609
610
611
    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);
612
613
614
615
616
617
618
619
620
    }
    CudaNonbondedMethod method = NO_CUTOFF;
    if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
        method = CUTOFF;
    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic) {
        method = PERIODIC;
    }
    data.customNonbondedMethod = method;

621
622
623
624
625
626
627
628
629
630
631
    // 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);
    }

632
633
634
    // Record information for the expressions.

    vector<string> paramNames;
635
636
    for (int i = 0; i < force.getNumPerParticleParameters(); i++)
        paramNames.push_back(force.getPerParticleParameterName(i));
637
638
639
640
    globalParamNames.resize(force.getNumGlobalParameters());
    globalParamValues.resize(force.getNumGlobalParameters());
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
        globalParamNames[i] = force.getGlobalParameterName(i);
641
        globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
642
    }
643
    gpuSetCustomNonbondedParameters(gpu, parameters, exclusionList, method, (float) force.getCutoffDistance(), force.getEnergyFunction(), paramNames, globalParamNames);
644
    if (globalParamValues.size() > 0)
Peter Eastman's avatar
Peter Eastman committed
645
        SetCustomNonbondedGlobalParams(globalParamValues);
646
647
648
}

void CudaCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
649
    updateGlobalParams(context);
650
651
652
}

double CudaCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
653
    updateGlobalParams(context);
654
655
656
    return 0.0;
}

657
658
void CudaCalcCustomNonbondedForceKernel::updateGlobalParams(ContextImpl& context) {
    bool changed = false;
659
    for (int i = 0; i < (int) globalParamNames.size(); i++) {
660
661
662
663
664
665
        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
666
        SetCustomNonbondedGlobalParams(globalParamValues);
667
668
}

669
CudaCalcGBSAOBCForceKernel::~CudaCalcGBSAOBCForceKernel() {
670
671
}

672
void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
673

Peter Eastman's avatar
Peter Eastman committed
674
    int numParticles = system.getNumParticles();
675
    _gpuContext* gpu = data.gpu;
Peter Eastman's avatar
Peter Eastman committed
676
677
    vector<float> radius(numParticles);
    vector<float> scale(numParticles);
678
    vector<float> charge(numParticles);
Peter Eastman's avatar
Peter Eastman committed
679
    for (int i = 0; i < numParticles; i++) {
680
681
        double particleCharge, particleRadius, scalingFactor;
        force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor);
Peter Eastman's avatar
Peter Eastman committed
682
        radius[i] = (float) particleRadius;
683
        scale[i] = (float) scalingFactor;
684
        charge[i] = (float) particleCharge;
685
    }
686
    gpuSetObcParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), radius, scale, charge);
687
688
}

689
void CudaCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
690
691
}

Mark Friedrichs's avatar
Mark Friedrichs committed
692
693
694
695
696
697
698
699
700
701
702
703
704
705
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++) {
706
        double charge, particleRadius, gamma;
Mark Friedrichs's avatar
Mark Friedrichs committed
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
        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;
}

724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
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
744
        SetCustomExternalGlobalParams(globalParamValues);
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
}

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
767
        SetCustomExternalGlobalParams(globalParamValues);
768
769
}

770
void OPENMMCUDA_EXPORT OpenMM::cudaOpenMMInitializeIntegration(const System& system, CudaPlatform::PlatformData& data, const Integrator& integrator) {
771
772
773
774
775
776
777
778
779
780
781
782
783
784

    // 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) {
785
786
        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>());
787
788
        if (gpu->bIncludeGBSA || gpu->bIncludeGBVI)
            throw OpenMMException("CudaPlatform requires GBSAOBCForce and GBVIForce to be used with a NonbondedForce");
789
    }
790
791
792
    
    // Set masses.
    
Peter Eastman's avatar
Peter Eastman committed
793
794
795
796
    int numParticles = system.getNumParticles();
    vector<float> mass(numParticles);
    for (int i = 0; i < numParticles; i++)
        mass[i] = (float) system.getParticleMass(i);
797
798
799
800
    gpuSetMass(gpu, mass);
    
    // Set constraints.
    
801
    int numConstraints = system.getNumConstraints();
Peter Eastman's avatar
Peter Eastman committed
802
803
    vector<int> particle1(numConstraints);
    vector<int> particle2(numConstraints);
804
805
806
807
    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
808
        int particle1Index, particle2Index;
809
        double constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
810
811
812
        system.getConstraintParameters(i, particle1Index, particle2Index, constraintDistance);
        particle1[i] = particle1Index;
        particle2[i] = particle2Index;
813
        distance[i] = (float) constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
814
815
        invMass1[i] = 1.0f/mass[particle1Index];
        invMass2[i] = 1.0f/mass[particle2Index];
816
    }
817
    gpuSetConstraintParameters(gpu, particle1, particle2, distance, invMass1, invMass2, (float)integrator.getConstraintTolerance());
818
819
    
    // Finish initialization.
820

821
822
823
    gpuBuildThreadBlockWorkList(gpu);
    gpuBuildExclusionList(gpu);
    gpuBuildOutputBuffers(gpu);
824
    gpuSetConstants(gpu);
825
826
827
828
    if (gpu->bIncludeGBSA || gpu->bIncludeGBVI)
        kClearBornSumAndForces(gpu);
    else
        kClearForces(gpu);
829
    cudaThreadSynchronize();
830
831
}

832
double CudaCalcGBSAOBCForceKernel::executeEnergy(ContextImpl& context) {
833
	return 0.0;
834
835
836
837
838
839
}

CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
}

void CudaIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
840
    cudaOpenMMInitializeIntegration(system, data, integrator);
841
842
843
    prevStepSize = -1.0;
}

844
void CudaIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
845
846
847
848
849
    _gpuContext* gpu = data.gpu;
    double stepSize = integrator.getStepSize();
    if (stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
850
        gpuSetVerletIntegrationParameters(gpu, (float) stepSize, 0.0f);
851
852
853
854
855
        gpuSetConstants(gpu);
        prevStepSize = stepSize;
    }
    kVerletUpdatePart1(gpu);
    kApplyFirstShake(gpu);
856
    kApplyFirstSettle(gpu);
857
    kApplyFirstCCMA(gpu);
858
859
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
860
861
            gpu->bCalculateCM = true;
    kVerletUpdatePart2(gpu);
862
    data.time += stepSize;
863
    data.stepCount++;
864
865
866
867
868
869
}

CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
}

void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
870
    cudaOpenMMInitializeIntegration(system, data, integrator);
871
872
873
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
874
875
    prevTemp = -1.0;
    prevFriction = -1.0;
876
877
878
    prevStepSize = -1.0;
}

879
void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
880
    _gpuContext* gpu = data.gpu;
881
882
883
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
884
885
    if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
886
        
887
        double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
888
        gpuSetLangevinIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature, 0.0f);
889
890
891
892
893
894
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
895
    kLangevinUpdatePart1(gpu);
896
897
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
898
            gpu->bCalculateCM = true;
899
    kLangevinUpdatePart2(gpu);
900
    kApplySecondShake(gpu);
901
    kApplySecondSettle(gpu);
902
    kApplySecondCCMA(gpu);
903
    kSetVelocitiesFromPositions(gpu);
904
    data.time += stepSize;
905
    data.stepCount++;
906
}
907
908
909
910
911

CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
}

void CudaIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
912
    cudaOpenMMInitializeIntegration(system, data, integrator);
913
914
915
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
916
917
    prevTemp = -1.0;
    prevFriction = -1.0;
918
919
920
    prevStepSize = -1.0;
}

921
void CudaIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
922
923
924
925
926
927
928
929
    _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);
930
        gpuSetBrownianIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature);
931
932
933
934
935
936
937
938
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
    kBrownianUpdatePart1(gpu);
    kApplyFirstShake(gpu);
939
    kApplyFirstSettle(gpu);
940
    kApplyFirstCCMA(gpu);
941
942
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
943
944
            gpu->bCalculateCM = true;
    kBrownianUpdatePart2(gpu);
945
    data.time += stepSize;
946
947
948
949
950
951
952
    data.stepCount++;
}

CudaIntegrateVariableVerletStepKernel::~CudaIntegrateVariableVerletStepKernel() {
}

void CudaIntegrateVariableVerletStepKernel::initialize(const System& system, const VariableVerletIntegrator& integrator) {
953
    cudaOpenMMInitializeIntegration(system, data, integrator);
954
955
956
    prevErrorTol = -1.0;
}

957
void CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
958
959
960
961
962
963
964
965
966
    _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;
    }
967
    float maxStepSize = (float)(maxTime-data.time);
968
    kSelectVerletStepSize(gpu, maxStepSize);
969
970
971
972
973
974
975
976
977
978
    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;
979
980
    if ((*gpu->psStepSize)[0].y == maxStepSize)
        data.time = maxTime; // Avoid round-off error
981
    data.stepCount++;
982
}
983

984
985
986
987
CudaIntegrateVariableLangevinStepKernel::~CudaIntegrateVariableLangevinStepKernel() {
}

void CudaIntegrateVariableLangevinStepKernel::initialize(const System& system, const VariableLangevinIntegrator& integrator) {
988
    cudaOpenMMInitializeIntegration(system, data, integrator);
989
990
991
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
992
993
    prevTemp = -1.0;
    prevFriction = -1.0;
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
    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);
1006
        gpuSetLangevinIntegrationParameters(gpu, (float) tau, 0.0f, (float) temperature, (float) errorTol);
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
        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);
1023
    kSetVelocitiesFromPositions(gpu);
1024
1025
1026
1027
1028
1029
1030
    gpu->psStepSize->Download();
    data.time += (*gpu->psStepSize)[0].y;
    if ((*gpu->psStepSize)[0].y == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
}

1031
1032
1033
1034
CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
}

void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
1035
1036
1037
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) thermostat.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
1038
1039
    prevTemp = -1.0;
    prevFrequency = -1.0;
1040
1041
1042
    prevStepSize = -1.0;
}

1043
void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
1044
    _gpuContext* gpu = data.gpu;
1045
1046
    double temperature = context.getParameter(AndersenThermostat::Temperature());
    double frequency = context.getParameter(AndersenThermostat::CollisionFrequency());
1047
1048
1049
1050
    double stepSize = context.getIntegrator().getStepSize();
    if (temperature != prevTemp || frequency != prevFrequency || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
1051
        gpuSetAndersenThermostatParameters(gpu, (float) temperature, (float) frequency);
1052
1053
1054
1055
1056
1057
1058
1059
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFrequency = frequency;
        prevStepSize = stepSize;
    }
    kCalculateAndersenThermostat(gpu);
}
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
1086
1087
1088
1089
1090
1091
1092
1093
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);
1094
1095
    for (int i = 0; i < (int) gpu->posCellOffsets.size(); i++)
        gpu->posCellOffsets[i] = make_int3(0, 0, 0);
1096
1097
1098
1099
1100
1101
1102
}

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

1103
void CudaCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
1104
1105
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
1106
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1107
        masses[i] = system.getParticleMass(i);
1108
1109
}

1110
double CudaCalcKineticEnergyKernel::execute(ContextImpl& context) {
1111
1112
1113
    // We don't currently have a GPU kernel to do this, so we retrieve the velocities and calculate the energy
    // on the CPU.
    
1114
1115
    _gpuContext* gpu = data.gpu;
    gpu->psVelm4->Download();
1116
    double energy = 0.0;
1117
1118
1119
1120
    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);
    }
1121
1122
    return 0.5*energy;
}
1123

1124
void CudaRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
1125
    data.removeCM = true;
1126
    data.cmMotionFrequency = force.getFrequency();
1127
1128
}

1129
void CudaRemoveCMMotionKernel::execute(ContextImpl& context) {
1130
}