CudaKernels.cpp 29.9 KB
Newer Older
1
2
3
4
5
6
7
8
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
9
 * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
10
11
12
 * Authors: Peter Eastman                                                     *
 * Contributors:                                                              *
 *                                                                            *
13
14
15
16
 * This program is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU Lesser General Public License as published   *
 * by the Free Software Foundation, either version 3 of the License, or       *
 * (at your option) any later version.                                        *
17
 *                                                                            *
18
19
20
21
 * This program is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
 * GNU Lesser General Public License for more details.                        *
22
 *                                                                            *
23
24
 * You should have received a copy of the GNU Lesser General Public License   *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
25
26
27
28
 * -------------------------------------------------------------------------- */

#include "CudaKernels.h"
#include "CudaStreamImpl.h"
29
#include "openmm/LangevinIntegrator.h"
30
#include "openmm/Context.h"
31
#include "ReferencePlatform.h"
32
#include "openmm/internal/ContextImpl.h"
33
#include "kernels/gputypes.h"
34
#include "kernels/cudaKernels.h"
35
36
37
38
39
40
41
#include <cmath>

extern "C" int gpuSetConstants( gpuContext gpu );

using namespace OpenMM;
using namespace std;

42
static void calcForces(ContextImpl& context, CudaPlatform::PlatformData& data) {
43
    _gpuContext* gpu = data.gpu;
44
    if (data.nonbondedMethod != NO_CUTOFF && data.computeForceCount%100 == 0)
45
        gpuReorderAtoms(gpu);
46
    data.computeForceCount++;
47
48
    kClearForces(gpu);
    if (gpu->bIncludeGBSA) {
49
        gpu->bRecalculateBornRadii = true;
50
51
52
53
        kCalculateCDLJObcGbsaForces1(gpu);
        kReduceObcGbsaBornForces(gpu);
        kCalculateObcGbsaForces2(gpu);
    }
54
    else if (data.hasNonbonded)
55
        kCalculateCDLJForces(gpu);
56
57
    if (data.hasCustomNonbonded)
        kCalculateCustomNonbondedForces(gpu, data.hasNonbonded);
58
    kCalculateLocalForces(gpu);
59
    kReduceForces(gpu);
60
61
}

62
static double calcEnergy(ContextImpl& context, System& system) {
63
64
    // We don't currently have GPU kernels to calculate energy, so instead we have the reference
    // platform do it.  This is VERY slow.
65

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

80
81
82
void CudaInitializeForcesKernel::initialize(const System& system) {
}

83
void CudaInitializeForcesKernel::execute(ContextImpl& context) {
84
85
}

86
87
88
void CudaUpdateTimeKernel::initialize(const System& system) {
}

89
double CudaUpdateTimeKernel::getTime(const ContextImpl& context) const {
90
91
92
    return data.time;
}

93
void CudaUpdateTimeKernel::setTime(ContextImpl& context, double time) {
94
95
96
    data.time = time;
}

97
98
99
100
101
102
103
CudaCalcHarmonicBondForceKernel::~CudaCalcHarmonicBondForceKernel() {
}

void CudaCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasBonds = true;
104
    numBonds = force.getNumBonds();
Peter Eastman's avatar
Peter Eastman committed
105
106
    vector<int> particle1(numBonds);
    vector<int> particle2(numBonds);
107
108
109
110
    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
111
        force.getBondParameters(i, particle1[i], particle2[i], lengthValue, kValue);
112
113
114
        length[i] = (float) lengthValue;
        k[i] = (float) kValue;
    }
Peter Eastman's avatar
Peter Eastman committed
115
    gpuSetBondParameters(data.gpu, particle1, particle2, length, k);
116
117
}

118
void CudaCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
119
120
121
122
    if (data.primaryKernel == this)
        calcForces(context, data);
}

123
double CudaCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
124
125
126
127
128
129
130
131
132
133
134
135
    if (data.primaryKernel == this)
        return calcEnergy(context, system);
    return 0.0;
}

CudaCalcHarmonicAngleForceKernel::~CudaCalcHarmonicAngleForceKernel() {
}

void CudaCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasAngles = true;
136
    numAngles = force.getNumAngles();
137
    const float RadiansToDegrees = (float) (180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
138
139
140
    vector<int> particle1(numAngles);
    vector<int> particle2(numAngles);
    vector<int> particle3(numAngles);
141
142
143
144
    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
145
        force.getAngleParameters(i, particle1[i], particle2[i], particle3[i], angleValue, kValue);
146
147
        angle[i] = (float) (angleValue*RadiansToDegrees);
        k[i] = (float) kValue;
148
    }
Peter Eastman's avatar
Peter Eastman committed
149
    gpuSetBondAngleParameters(data.gpu, particle1, particle2, particle3, angle, k);
150
}
151

152
void CudaCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
153
154
155
156
    if (data.primaryKernel == this)
        calcForces(context, data);
}

157
double CudaCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
158
159
160
161
162
163
164
165
166
167
168
169
170
    if (data.primaryKernel == this)
        return calcEnergy(context, system);
    return 0.0;
}

CudaCalcPeriodicTorsionForceKernel::~CudaCalcPeriodicTorsionForceKernel() {
}

void CudaCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasPeriodicTorsions = true;
    numTorsions = force.getNumTorsions();
171
    const float RadiansToDegrees = (float)(180.0/3.14159265);
Peter Eastman's avatar
Peter Eastman committed
172
173
174
175
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
176
177
178
179
180
    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
181
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], periodicity[i], phaseValue, kValue);
182
183
        k[i] = (float) kValue;
        phase[i] = (float) (phaseValue*RadiansToDegrees);
184
    }
Peter Eastman's avatar
Peter Eastman committed
185
    gpuSetDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, k, phase, periodicity);
186
187
}

188
void CudaCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
189
190
191
192
    if (data.primaryKernel == this)
        calcForces(context, data);
}

193
double CudaCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
194
195
196
197
198
199
200
201
202
203
204
205
206
    if (data.primaryKernel == this)
        return calcEnergy(context, system);
    return 0.0;
}

CudaCalcRBTorsionForceKernel::~CudaCalcRBTorsionForceKernel() {
}

void CudaCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasRB = true;
    numTorsions = force.getNumTorsions();
Peter Eastman's avatar
Peter Eastman committed
207
208
209
210
    vector<int> particle1(numTorsions);
    vector<int> particle2(numTorsions);
    vector<int> particle3(numTorsions);
    vector<int> particle4(numTorsions);
211
212
213
214
215
216
217
218
    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
219
        force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], c[0], c[1], c[2], c[3], c[4], c[5]);
220
221
222
223
224
225
        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];
226
    }
Peter Eastman's avatar
Peter Eastman committed
227
    gpuSetRbDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
228
229
}

230
void CudaCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
231
232
233
234
    if (data.primaryKernel == this)
        calcForces(context, data);
}

235
double CudaCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
236
237
238
239
240
241
242
243
    if (data.primaryKernel == this)
        return calcEnergy(context, system);
    return 0.0;
}

CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() {
}

244
void CudaCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
245
246
247
    if (data.primaryKernel == NULL)
        data.primaryKernel = this;
    data.hasNonbonded = true;
Peter Eastman's avatar
Peter Eastman committed
248
    numParticles = force.getNumParticles();
249
    _gpuContext* gpu = data.gpu;
250
251
252
253

    // Identify which exceptions are 1-4 interactions.

    vector<pair<int, int> > exclusions;
254
    vector<int> exceptions;
255
256
257
258
259
260
    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)
261
            exceptions.push_back(i);
262
263
    }

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

        if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
            Vec3 boxVectors[3];
            force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
301
            gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
302
303
304
305
306
307
            double ewaldErrorTol = force.getEwaldErrorTolerance();
            double alpha = (1.0/force.getCutoffDistance())*std::sqrt(-std::log(ewaldErrorTol));
            double mx = boxVectors[0][0]/force.getCutoffDistance();
            double my = boxVectors[1][1]/force.getCutoffDistance();
            double mz = boxVectors[2][2]/force.getCutoffDistance();
            double pi = 3.1415926535897932385;
308
309
310
            int kmaxx = (int)std::ceil(-(mx/pi)*std::log(ewaldErrorTol));
            int kmaxy = (int)std::ceil(-(my/pi)*std::log(ewaldErrorTol));
            int kmaxz = (int)std::ceil(-(mz/pi)*std::log(ewaldErrorTol));
311
312
313
314
315
316
            if (kmaxx%2 == 0)
                kmaxx++;
            if (kmaxy%2 == 0)
                kmaxy++;
            if (kmaxz%2 == 0)
                kmaxz++;
317
            gpuSetEwaldParameters(gpu, (float)alpha, kmaxx, kmaxy, kmaxz);
318
319
            method = EWALD;
        }
320
321
        data.nonbondedMethod = method;
        gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method);
322
323
324
325
    }

    // Initialize 1-4 nonbonded interactions.
    
326
    {
327
328
329
330
331
332
333
334
        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++) {
335
            double charge, sig, eps;
336
            force.getExceptionParameters(exceptions[i], particle1[i], particle2[i], charge, sig, eps);
337
338
            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
339
340
            q1[i] = (float) charge;
            q2[i] = 1.0f;
341
        }
Peter Eastman's avatar
Peter Eastman committed
342
        gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, particle1, particle2, c6, c12, q1, q2);
343
344
345
    }
}

346
void CudaCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
347
348
    if (data.primaryKernel == this)
        calcForces(context, data);
349
350
}

351
double CudaCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
352
353
354
    if (data.primaryKernel == this)
        return calcEnergy(context, system);
    return 0.0;
355
356
}

357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
CudaCalcCustomNonbondedForceKernel::~CudaCalcCustomNonbondedForceKernel() {
}

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

    // Identify which exceptions are actual interactions.

    vector<pair<int, int> > exclusions;
    vector<int> exceptions;
    {
        vector<double> parameters;
        for (int i = 0; i < force.getNumExceptions(); i++) {
            int particle1, particle2;
            force.getExceptionParameters(i, particle1, particle2, parameters);
            exclusions.push_back(pair<int, int>(particle1, particle2));
            if (parameters.size() > 0)
                exceptions.push_back(i);
        }
    }

    // Initialize nonbonded interactions.

    vector<int> particle(numParticles);
    vector<vector<double> > parameters(numParticles);
    vector<vector<int> > exclusionList(numParticles);
    for (int i = 0; i < numParticles; i++) {
        force.getParticleParameters(i, parameters[i]);
        particle[i] = i;
        exclusionList[i].push_back(i);
    }
    for (int i = 0; i < (int)exclusions.size(); i++) {
        exclusionList[exclusions[i].first].push_back(exclusions[i].second);
        exclusionList[exclusions[i].second].push_back(exclusions[i].first);
    }
    CudaNonbondedMethod method = NO_CUTOFF;
    if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
        method = CUTOFF;
    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic) {
        Vec3 boxVectors[3];
        force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
        gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
        method = PERIODIC;
    }
    data.customNonbondedMethod = method;

    // Initialize exceptions.

    int numExceptions = exceptions.size();
    vector<int> exceptionParticle1(numExceptions);
    vector<int> exceptionParticle2(numExceptions);
    vector<vector<double> > exceptionParams(numExceptions);
    for (int i = 0; i < numExceptions; i++)
        force.getExceptionParameters(exceptions[i], exceptionParticle1[i], exceptionParticle2[i], exceptionParams[i]);

    // Record information for the expressions.

    vector<string> paramNames;
    vector<string> combiningRules;
    for (int i = 0; i < force.getNumParameters(); i++) {
        paramNames.push_back(force.getParameterName(i));
        combiningRules.push_back(force.getParameterCombiningRule(i));
    }
    gpuSetCustomNonbondedParameters(gpu, parameters, exclusionList, exceptionParticle1, exceptionParticle2, exceptionParams, method,
            (float)force.getCutoffDistance(), force.getEnergyFunction(), combiningRules, paramNames);
}

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

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

439
CudaCalcGBSAOBCForceKernel::~CudaCalcGBSAOBCForceKernel() {
440
441
}

442
void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
443

Peter Eastman's avatar
Peter Eastman committed
444
    int numParticles = system.getNumParticles();
445
    _gpuContext* gpu = data.gpu;
Peter Eastman's avatar
Peter Eastman committed
446
447
    vector<float> radius(numParticles);
    vector<float> scale(numParticles);
448
    vector<float> charge(numParticles);
Peter Eastman's avatar
Peter Eastman committed
449
    for (int i = 0; i < numParticles; i++) {
450
451
        double particleCharge, particleRadius, scalingFactor;
        force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor);
Peter Eastman's avatar
Peter Eastman committed
452
        radius[i] = (float) particleRadius;
453
        scale[i] = (float) scalingFactor;
454
        charge[i] = (float) particleCharge;
455
    }
456
    gpuSetObcParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), radius, scale, charge);
457
458
}

459
void CudaCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
460
461
}

462
static void initializeIntegration(const System& system, CudaPlatform::PlatformData& data, const Integrator& integrator) {
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479

    // Initialize any terms that haven't already been handled by a Force.

    _gpuContext* gpu = data.gpu;
    if (!data.hasBonds)
        gpuSetBondParameters(gpu, vector<int>(), vector<int>(), vector<float>(), vector<float>());
    if (!data.hasAngles)
        gpuSetBondAngleParameters(gpu, vector<int>(), vector<int>(), vector<int>(), vector<float>(), vector<float>());
    if (!data.hasPeriodicTorsions)
        gpuSetDihedralParameters(gpu, vector<int>(), vector<int>(), vector<int>(), vector<int>(), vector<float>(), vector<float>(), vector<int>());
    if (!data.hasRB)
        gpuSetRbDihedralParameters(gpu, vector<int>(), vector<int>(), vector<int>(), vector<int>(), vector<float>(), vector<float>(),
                vector<float>(), vector<float>(), vector<float>(), vector<float>());
    if (!data.hasNonbonded) {
        gpuSetCoulombParameters(gpu, 138.935485f, vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<char>(), vector<vector<int> >(), NO_CUTOFF);
        gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, vector<int>(), vector<int>(), vector<float>(), vector<float>(), vector<float>(), vector<float>());
    }
480
481
482
    
    // Set masses.
    
Peter Eastman's avatar
Peter Eastman committed
483
484
485
486
    int numParticles = system.getNumParticles();
    vector<float> mass(numParticles);
    for (int i = 0; i < numParticles; i++)
        mass[i] = (float) system.getParticleMass(i);
487
488
489
490
    gpuSetMass(gpu, mass);
    
    // Set constraints.
    
491
    int numConstraints = system.getNumConstraints();
Peter Eastman's avatar
Peter Eastman committed
492
493
    vector<int> particle1(numConstraints);
    vector<int> particle2(numConstraints);
494
495
496
497
    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
498
        int particle1Index, particle2Index;
499
        double constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
500
501
502
        system.getConstraintParameters(i, particle1Index, particle2Index, constraintDistance);
        particle1[i] = particle1Index;
        particle2[i] = particle2Index;
503
        distance[i] = (float) constraintDistance;
Peter Eastman's avatar
Peter Eastman committed
504
505
        invMass1[i] = 1.0f/mass[particle1Index];
        invMass2[i] = 1.0f/mass[particle2Index];
506
    }
507
    gpuSetConstraintParameters(gpu, particle1, particle2, distance, invMass1, invMass2, (float)integrator.getConstraintTolerance());
508
509
    
    // Finish initialization.
510

511
512
513
    gpuBuildThreadBlockWorkList(gpu);
    gpuBuildExclusionList(gpu);
    gpuBuildOutputBuffers(gpu);
514
    gpuSetConstants(gpu);
515
516
517
    kClearBornForces(gpu);
    kClearForces(gpu);
    cudaThreadSynchronize();
518
519
}

520
double CudaCalcGBSAOBCForceKernel::executeEnergy(ContextImpl& context) {
521
	return 0.0;
522
523
524
525
526
527
}

CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
}

void CudaIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
528
    initializeIntegration(system, data, integrator);
529
530
531
    prevStepSize = -1.0;
}

532
void CudaIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
533
534
535
536
537
    _gpuContext* gpu = data.gpu;
    double stepSize = integrator.getStepSize();
    if (stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
538
        gpuSetVerletIntegrationParameters(gpu, (float) stepSize, 0.0f);
539
540
541
542
543
        gpuSetConstants(gpu);
        prevStepSize = stepSize;
    }
    kVerletUpdatePart1(gpu);
    kApplyFirstShake(gpu);
544
    kApplyFirstSettle(gpu);
545
    kApplyFirstCCMA(gpu);
546
547
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
548
549
            gpu->bCalculateCM = true;
    kVerletUpdatePart2(gpu);
550
    data.time += stepSize;
551
    data.stepCount++;
552
553
554
555
556
557
}

CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
}

void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
558
    initializeIntegration(system, data, integrator);
559
560
561
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
562
563
564
    prevStepSize = -1.0;
}

565
void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
566
    _gpuContext* gpu = data.gpu;
567
568
569
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
570
571
572
573
    if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
        double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
574
        gpuSetLangevinIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature, 0.0f);
575
576
577
578
579
580
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
581
    kLangevinUpdatePart1(gpu);
582
    kApplyFirstShake(gpu);
583
    kApplyFirstSettle(gpu);
584
    kApplyFirstCCMA(gpu);
585
586
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
587
            gpu->bCalculateCM = true;
588
    kLangevinUpdatePart2(gpu);
589
    kApplySecondShake(gpu);
590
    kApplySecondSettle(gpu);
591
    kApplySecondCCMA(gpu);
592
    data.time += stepSize;
593
    data.stepCount++;
594
}
595
596
597
598
599

CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
}

void CudaIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
600
    initializeIntegration(system, data, integrator);
601
602
603
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
604
605
606
    prevStepSize = -1.0;
}

607
void CudaIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
608
609
610
611
612
613
614
615
    _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);
616
        gpuSetBrownianIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature);
617
618
619
620
621
622
623
624
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
    kBrownianUpdatePart1(gpu);
    kApplyFirstShake(gpu);
625
    kApplyFirstSettle(gpu);
626
    kApplyFirstCCMA(gpu);
627
628
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
629
630
            gpu->bCalculateCM = true;
    kBrownianUpdatePart2(gpu);
631
    data.time += stepSize;
632
633
634
635
636
637
638
639
640
641
642
    data.stepCount++;
}

CudaIntegrateVariableVerletStepKernel::~CudaIntegrateVariableVerletStepKernel() {
}

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

643
void CudaIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
644
645
646
647
648
649
650
651
652
    _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;
    }
653
    float maxStepSize = (float)(maxTime-data.time);
654
    kSelectVerletStepSize(gpu, maxStepSize);
655
656
657
658
659
660
661
662
663
664
    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;
665
666
    if ((*gpu->psStepSize)[0].y == maxStepSize)
        data.time = maxTime; // Avoid round-off error
667
    data.stepCount++;
668
}
669

670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
CudaIntegrateVariableLangevinStepKernel::~CudaIntegrateVariableLangevinStepKernel() {
}

void CudaIntegrateVariableLangevinStepKernel::initialize(const System& system, const VariableLangevinIntegrator& integrator) {
    initializeIntegration(system, data, integrator);
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) integrator.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
    prevErrorTol = -1.0;
}

void CudaIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
    _gpuContext* gpu = data.gpu;
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double errorTol = integrator.getErrorTolerance();
    if (temperature != prevTemp || friction != prevFriction || errorTol != prevErrorTol) {
        // Initialize the GPU parameters.

        double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
        gpuSetLangevinIntegrationParameters(gpu, (float) tau, 0.0f, (float) temperature, errorTol);
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFriction = friction;
        prevErrorTol = errorTol;
    }
    float maxStepSize = (float)(maxTime-data.time);
    kSelectLangevinStepSize(gpu, maxStepSize);
    kLangevinUpdatePart1(gpu);
    kApplyFirstShake(gpu);
    kApplyFirstSettle(gpu);
    kApplyFirstCCMA(gpu);
    if (data.removeCM)
        if (data.stepCount%data.cmMotionFrequency == 0)
            gpu->bCalculateCM = true;
    kLangevinUpdatePart2(gpu);
    kApplySecondShake(gpu);
    kApplySecondSettle(gpu);
    kApplySecondCCMA(gpu);
    gpu->psStepSize->Download();
    data.time += (*gpu->psStepSize)[0].y;
    if ((*gpu->psStepSize)[0].y == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
}

717
718
719
720
CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
}

void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
721
722
723
    _gpuContext* gpu = data.gpu;
    gpu->seed = (unsigned long) thermostat.getRandomNumberSeed();
    gpuInitializeRandoms(gpu);
724
725
726
    prevStepSize = -1.0;
}

727
void CudaApplyAndersenThermostatKernel::execute(ContextImpl& context) {
728
    _gpuContext* gpu = data.gpu;
729
730
    double temperature = context.getParameter(AndersenThermostat::Temperature());
    double frequency = context.getParameter(AndersenThermostat::CollisionFrequency());
731
732
733
734
    double stepSize = context.getIntegrator().getStepSize();
    if (temperature != prevTemp || frequency != prevFrequency || stepSize != prevStepSize) {
        // Initialize the GPU parameters.
        
735
        gpuSetAndersenThermostatParameters(gpu, (float) temperature, frequency);
736
737
738
739
740
741
742
743
        gpuSetConstants(gpu);
        kGenerateRandoms(gpu);
        prevTemp = temperature;
        prevFrequency = frequency;
        prevStepSize = stepSize;
    }
    kCalculateAndersenThermostat(gpu);
}
744

745
void CudaCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
746
747
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
748
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
749
        masses[i] = system.getParticleMass(i);
750
751
}

752
double CudaCalcKineticEnergyKernel::execute(ContextImpl& context) {
753
754
755
    // We don't currently have a GPU kernel to do this, so we retrieve the velocities and calculate the energy
    // on the CPU.
    
756
    const Stream& velocities = context.getVelocities();
757
758
759
760
761
762
763
764
    double* v = new double[velocities.getSize()*3];
    velocities.saveToArray(v);
    double energy = 0.0;
    for (size_t i = 0; i < masses.size(); ++i)
        energy += masses[i]*(v[i*3]*v[i*3]+v[i*3+1]*v[i*3+1]+v[i*3+2]*v[i*3+2]);
    delete v;
    return 0.5*energy;
}
765

766
void CudaRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
767
    data.removeCM = true;
768
    data.cmMotionFrequency = force.getFrequency();
769
770
}

771
void CudaRemoveCMMotionKernel::execute(ContextImpl& context) {
772
}