CudaKernels.cpp 43.6 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
CudaCalcHarmonicBondForceKernel::~CudaCalcHarmonicBondForceKernel() {
}

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

221
void CudaCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
222
223
}

224
double CudaCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
225
226
227
    return 0.0;
}

228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
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
249
        SetCustomBondGlobalParams(globalParamValues);
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
}

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
272
        SetCustomBondGlobalParams(globalParamValues);
273
274
}

275
276
277
278
279
CudaCalcHarmonicAngleForceKernel::~CudaCalcHarmonicAngleForceKernel() {
}

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

296
void CudaCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
297
298
}

299
double CudaCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
300
301
302
    return 0.0;
}

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
348
349
350
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);
}

351
352
353
354
355
356
CudaCalcPeriodicTorsionForceKernel::~CudaCalcPeriodicTorsionForceKernel() {
}

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

374
void CudaCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
375
376
}

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

410
void CudaCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
411
412
}

413
double CudaCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
414
415
416
    return 0.0;
}

417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
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);
}

466
467
468
CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() {
}

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

    // Identify which exceptions are 1-4 interactions.

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

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

        // Compute the Ewald self energy.

        data.ewaldSelfEnergy = 0.0;
539
540
541
542
543
        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];
        }
544
545
546
547
548
549
550

        // Compute the long range dispersion correction.

        if (force.getUseDispersionCorrection())
            data.dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
        else
            data.dispersionCoefficient = 0.0;
551
552
553
554
    }

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

575
void CudaCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
576
577
}

578
double CudaCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
579
    return 0.0;
580
581
}

582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
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);
    }
600
601
602
603
604
    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);
605
606
607
608
609
610
611
612
613
    }
    CudaNonbondedMethod method = NO_CUTOFF;
    if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
        method = CUTOFF;
    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic) {
        method = PERIODIC;
    }
    data.customNonbondedMethod = method;

614
615
616
617
618
619
620
621
622
623
624
    // 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);
    }

625
626
627
    // Record information for the expressions.

    vector<string> paramNames;
628
629
    for (int i = 0; i < force.getNumPerParticleParameters(); i++)
        paramNames.push_back(force.getPerParticleParameterName(i));
630
631
632
633
    globalParamNames.resize(force.getNumGlobalParameters());
    globalParamValues.resize(force.getNumGlobalParameters());
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
        globalParamNames[i] = force.getGlobalParameterName(i);
634
        globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
635
    }
636
    gpuSetCustomNonbondedParameters(gpu, parameters, exclusionList, method, (float) force.getCutoffDistance(), force.getEnergyFunction(), paramNames, globalParamNames);
637
    if (globalParamValues.size() > 0)
Peter Eastman's avatar
Peter Eastman committed
638
        SetCustomNonbondedGlobalParams(globalParamValues);
639
640
641
}

void CudaCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
642
    updateGlobalParams(context);
643
644
645
}

double CudaCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
646
    updateGlobalParams(context);
647
648
649
    return 0.0;
}

650
651
void CudaCalcCustomNonbondedForceKernel::updateGlobalParams(ContextImpl& context) {
    bool changed = false;
652
    for (int i = 0; i < (int) globalParamNames.size(); i++) {
653
654
655
656
657
658
        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
659
        SetCustomNonbondedGlobalParams(globalParamValues);
660
661
}

662
CudaCalcGBSAOBCForceKernel::~CudaCalcGBSAOBCForceKernel() {
663
664
}

665
void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
666

Peter Eastman's avatar
Peter Eastman committed
667
    int numParticles = system.getNumParticles();
668
    _gpuContext* gpu = data.gpu;
Peter Eastman's avatar
Peter Eastman committed
669
670
    vector<float> radius(numParticles);
    vector<float> scale(numParticles);
671
    vector<float> charge(numParticles);
Peter Eastman's avatar
Peter Eastman committed
672
    for (int i = 0; i < numParticles; i++) {
673
674
        double particleCharge, particleRadius, scalingFactor;
        force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor);
Peter Eastman's avatar
Peter Eastman committed
675
        radius[i] = (float) particleRadius;
676
        scale[i] = (float) scalingFactor;
677
        charge[i] = (float) particleCharge;
678
    }
679
    gpuSetObcParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), radius, scale, charge);
680
681
}

682
void CudaCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
683
684
}

Mark Friedrichs's avatar
Mark Friedrichs committed
685
686
687
688
689
690
691
692
693
694
695
696
697
698
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++) {
699
        double charge, particleRadius, gamma;
Mark Friedrichs's avatar
Mark Friedrichs committed
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
        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;
}

717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
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
737
        SetCustomExternalGlobalParams(globalParamValues);
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
}

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
760
        SetCustomExternalGlobalParams(globalParamValues);
761
762
}

763
void OPENMMCUDA_EXPORT OpenMM::cudaOpenMMInitializeIntegration(const System& system, CudaPlatform::PlatformData& data, const Integrator& integrator) {
764
765
766
767
768
769
770
771
772
773
774
775
776
777

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

814
815
816
    gpuBuildThreadBlockWorkList(gpu);
    gpuBuildExclusionList(gpu);
    gpuBuildOutputBuffers(gpu);
817
    gpuSetConstants(gpu);
818
819
820
821
    if (gpu->bIncludeGBSA || gpu->bIncludeGBVI)
        kClearBornSumAndForces(gpu);
    else
        kClearForces(gpu);
822
    cudaThreadSynchronize();
823
824
}

825
double CudaCalcGBSAOBCForceKernel::executeEnergy(ContextImpl& context) {
826
	return 0.0;
827
828
829
830
831
832
}

CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
}

void CudaIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
833
    cudaOpenMMInitializeIntegration(system, data, integrator);
834
835
836
    prevStepSize = -1.0;
}

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

CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
}

void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
863
    cudaOpenMMInitializeIntegration(system, data, integrator);
864
865
866
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
867
868
    prevTemp = -1.0;
    prevFriction = -1.0;
869
870
871
    prevStepSize = -1.0;
}

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

CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
}

void CudaIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
905
    cudaOpenMMInitializeIntegration(system, data, integrator);
906
907
908
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
909
910
    prevTemp = -1.0;
    prevFriction = -1.0;
911
912
913
    prevStepSize = -1.0;
}

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

CudaIntegrateVariableVerletStepKernel::~CudaIntegrateVariableVerletStepKernel() {
}

void CudaIntegrateVariableVerletStepKernel::initialize(const System& system, const VariableVerletIntegrator& integrator) {
946
    cudaOpenMMInitializeIntegration(system, data, integrator);
947
948
949
    prevErrorTol = -1.0;
}

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

977
978
979
980
CudaIntegrateVariableLangevinStepKernel::~CudaIntegrateVariableLangevinStepKernel() {
}

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

1024
1025
1026
1027
CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
}

void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
1028
1029
1030
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) thermostat.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
1031
1032
    prevTemp = -1.0;
    prevFrequency = -1.0;
1033
1034
1035
    prevStepSize = -1.0;
}

1036
void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
1037
    _gpuContext* gpu = data.gpu;
1038
1039
    double temperature = context.getParameter(AndersenThermostat::Temperature());
    double frequency = context.getParameter(AndersenThermostat::CollisionFrequency());
1040
1041
1042
1043
    double stepSize = context.getIntegrator().getStepSize();
    if (temperature != prevTemp || frequency != prevFrequency || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
1044
        gpuSetAndersenThermostatParameters(gpu, (float) temperature, (float) frequency);
1045
1046
1047
1048
1049
1050
1051
1052
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFrequency = frequency;
        prevStepSize = stepSize;
    }
    kCalculateAndersenThermostat(gpu);
}
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
1086
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);
1087
1088
    for (int i = 0; i < (int) gpu->posCellOffsets.size(); i++)
        gpu->posCellOffsets[i] = make_int3(0, 0, 0);
1089
1090
1091
1092
1093
1094
1095
}

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

1096
void CudaCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
1097
1098
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
1099
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1100
        masses[i] = system.getParticleMass(i);
1101
1102
}

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

1117
void CudaRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
1118
    data.removeCM = true;
1119
    data.cmMotionFrequency = force.getFrequency();
1120
1121
}

1122
void CudaRemoveCMMotionKernel::execute(ContextImpl& context) {
1123
}