OpenCLKernels.h 58.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#ifndef OPENMM_OPENCLKERNELS_H_
#define OPENMM_OPENCLKERNELS_H_

/* -------------------------------------------------------------------------- *
 *                                   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.               *
 *                                                                            *
12
 * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
 * Authors: Peter Eastman                                                     *
 * Contributors:                                                              *
 *                                                                            *
 * 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.                                        *
 *                                                                            *
 * 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.                        *
 *                                                                            *
 * 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/>.      *
 * -------------------------------------------------------------------------- */

#include "OpenCLPlatform.h"
31
32
#include "OpenCLArray.h"
#include "OpenCLContext.h"
33
#include "OpenCLFFT3D.h"
34
#include "OpenCLParameterSet.h"
35
#include "OpenCLSort.h"
36
37
38
39
40
#include "openmm/kernels.h"
#include "openmm/System.h"

namespace OpenMM {

41
/**
42
43
44
 * This kernel is invoked at the beginning and end of force and energy computations.  It gives the
 * Platform a chance to clear buffers and do other initialization at the beginning, and to do any
 * necessary work at the end to determine the final results.
45
 */
46
class OpenCLCalcForcesAndEnergyKernel : public CalcForcesAndEnergyKernel {
47
public:
48
    OpenCLCalcForcesAndEnergyKernel(std::string name, const Platform& platform, OpenCLContext& cl) : CalcForcesAndEnergyKernel(name, platform), cl(cl) {
49
50
51
52
53
54
55
56
    }
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     */
    void initialize(const System& system);
    /**
57
     * This is called at the beginning of each force/energy computation, before calcForcesAndEnergy() has been called on
58
     * any ForceImpl.
59
     *
60
61
62
     * @param context       the context in which to execute this kernel
     * @param includeForce  true if forces should be computed
     * @param includeEnergy true if potential energy should be computed
63
     * @param groups        a set of bit flags for which force groups to include
64
     */
65
    void beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
66
    /**
67
     * This is called at the end of each force/energy computation, after calcForcesAndEnergy() has been called on
68
69
     * every ForceImpl.
     *
70
71
72
     * @param context       the context in which to execute this kernel
     * @param includeForce  true if forces should be computed
     * @param includeEnergy true if potential energy should be computed
73
     * @param groups        a set of bit flags for which force groups to include
74
75
     * @param valid         the method may set this to false to indicate the results are invalid and the force/energy
     *                      calculation should be repeated
76
     * @return the potential energy of the system.  This value is added to all values returned by ForceImpls'
77
     * calcForcesAndEnergy() methods.  That is, each force kernel may <i>either</i> return its contribution to the
78
79
     * energy directly, <i>or</i> add it to an internal buffer so that it will be included here.
     */
80
    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups, bool& valid);
81
private:
82
   OpenCLContext& cl;
83
84
85
};

/**
86
87
 * This kernel provides methods for setting and retrieving various state data: time, positions,
 * velocities, and forces.
88
 */
89
class OpenCLUpdateStateDataKernel : public UpdateStateDataKernel {
90
public:
91
    OpenCLUpdateStateDataKernel(std::string name, const Platform& platform, OpenCLContext& cl) : UpdateStateDataKernel(name, platform), cl(cl) {
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
    }
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     */
    void initialize(const System& system);
    /**
     * Get the current time (in picoseconds).
     *
     * @param context    the context in which to execute this kernel
     */
    double getTime(const ContextImpl& context) const;
    /**
     * Set the current time (in picoseconds).
     *
     * @param context    the context in which to execute this kernel
     */
    void setTime(ContextImpl& context, double time);
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
    /**
     * Get the positions of all particles.
     *
     * @param positions  on exit, this contains the particle positions
     */
    void getPositions(ContextImpl& context, std::vector<Vec3>& positions);
    /**
     * Set the positions of all particles.
     *
     * @param positions  a vector containg the particle positions
     */
    void setPositions(ContextImpl& context, const std::vector<Vec3>& positions);
    /**
     * Get the velocities of all particles.
     *
     * @param velocities  on exit, this contains the particle velocities
     */
    void getVelocities(ContextImpl& context, std::vector<Vec3>& velocities);
    /**
     * Set the velocities of all particles.
     *
     * @param velocities  a vector containg the particle velocities
     */
    void setVelocities(ContextImpl& context, const std::vector<Vec3>& velocities);
    /**
     * Get the current forces on all particles.
     *
     * @param forces  on exit, this contains the forces
     */
    void getForces(ContextImpl& context, std::vector<Vec3>& forces);
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
    /**
     * Get the current periodic box vectors.
     *
     * @param a      on exit, this contains the vector defining the first edge of the periodic box
     * @param b      on exit, this contains the vector defining the second edge of the periodic box
     * @param c      on exit, this contains the vector defining the third edge of the periodic box
     */
    void getPeriodicBoxVectors(ContextImpl& context, Vec3& a, Vec3& b, Vec3& c) const;
    /**
     * Set the current periodic box vectors.
     *
     * @param a      the vector defining the first edge of the periodic box
     * @param b      the vector defining the second edge of the periodic box
     * @param c      the vector defining the third edge of the periodic box
     */
    void setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c) const;
Peter Eastman's avatar
Peter Eastman committed
157
158
159
160
161
162
163
164
165
166
167
168
    /**
     * Create a checkpoint recording the current state of the Context.
     * 
     * @param stream    an output stream the checkpoint data should be written to
     */
    void createCheckpoint(ContextImpl& context, std::ostream& stream);
    /**
     * Load a checkpoint that was written by createCheckpoint().
     * 
     * @param stream    an input stream the checkpoint data should be read from
     */
    void loadCheckpoint(ContextImpl& context, std::istream& stream);
169
private:
170
    OpenCLContext& cl;
171
};
172

173
174
175
176
177
/**
 * This kernel modifies the positions of particles to enforce distance constraints.
 */
class OpenCLApplyConstraintsKernel : public ApplyConstraintsKernel {
public:
178
179
    OpenCLApplyConstraintsKernel(std::string name, const Platform& platform, OpenCLContext& cl) : ApplyConstraintsKernel(name, platform),
            cl(cl), hasInitializedKernel(false) {
180
181
182
183
184
185
186
187
188
189
190
191
192
193
    }
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     */
    void initialize(const System& system);
    /**
     * Update particle positions to enforce constraints.
     *
     * @param context    the context in which to execute this kernel
     * @param tol        the distance tolerance within which constraints must be satisfied.
     */
    void apply(ContextImpl& context, double tol);
194
195
196
197
198
199
200
    /**
     * Update particle velocities to enforce constraints.
     *
     * @param context    the context in which to execute this kernel
     * @param tol        the velocity tolerance within which constraints must be satisfied.
     */
    void applyToVelocities(ContextImpl& context, double tol);
201
202
private:
    OpenCLContext& cl;
203
204
    bool hasInitializedKernel;
    cl::Kernel applyDeltasKernel;
205
206
};

207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
/**
 * This kernel recomputes the positions of virtual sites.
 */
class OpenCLVirtualSitesKernel : public VirtualSitesKernel {
public:
    OpenCLVirtualSitesKernel(std::string name, const Platform& platform, OpenCLContext& cl) : VirtualSitesKernel(name, platform), cl(cl) {
    }
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     */
    void initialize(const System& system);
    /**
     * Compute the virtual site locations.
     *
     * @param context    the context in which to execute this kernel
     */
    void computePositions(ContextImpl& context);
private:
    OpenCLContext& cl;
};

230
231
232
233
234
/**
 * This kernel is invoked by HarmonicBondForce to calculate the forces acting on the system and the energy of the system.
 */
class OpenCLCalcHarmonicBondForceKernel : public CalcHarmonicBondForceKernel {
public:
235
    OpenCLCalcHarmonicBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcHarmonicBondForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
236
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
237
238
239
240
241
242
243
244
245
246
    }
    ~OpenCLCalcHarmonicBondForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the HarmonicBondForce this kernel will be used for
     */
    void initialize(const System& system, const HarmonicBondForce& force);
    /**
247
     * Execute the kernel to calculate the forces and/or energy.
248
     *
249
250
251
252
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
253
     */
254
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
255
256
257
258
259
260
261
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the HarmonicBondForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const HarmonicBondForce& force);
262
263
private:
    int numBonds;
264
    bool hasInitializedKernel;
265
    OpenCLContext& cl;
266
    const System& system;
267
    OpenCLArray* params;
268
269
};

270
271
272
273
274
/**
 * This kernel is invoked by CustomBondForce to calculate the forces acting on the system and the energy of the system.
 */
class OpenCLCalcCustomBondForceKernel : public CalcCustomBondForceKernel {
public:
275
    OpenCLCalcCustomBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomBondForceKernel(name, platform),
276
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
277
278
279
280
281
282
283
284
285
286
    }
    ~OpenCLCalcCustomBondForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomBondForce this kernel will be used for
     */
    void initialize(const System& system, const CustomBondForce& force);
    /**
287
     * Execute the kernel to calculate the forces and/or energy.
288
     *
289
290
291
292
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
293
     */
294
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
295
296
297
298
299
300
301
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomBondForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomBondForce& force);
302
303
304
305
private:
    int numBonds;
    bool hasInitializedKernel;
    OpenCLContext& cl;
306
    const System& system;
307
    OpenCLParameterSet* params;
308
    OpenCLArray* globals;
309
310
311
312
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

313
314
315
316
317
/**
 * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
 */
class OpenCLCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
public:
318
    OpenCLCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcHarmonicAngleForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
319
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
320
321
322
323
324
325
326
327
328
329
    }
    ~OpenCLCalcHarmonicAngleForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the HarmonicAngleForce this kernel will be used for
     */
    void initialize(const System& system, const HarmonicAngleForce& force);
    /**
330
     * Execute the kernel to calculate the forces and/or energy.
331
     *
332
333
334
335
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
336
     */
337
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
338
339
340
341
342
343
344
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the HarmonicAngleForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force);
345
346
private:
    int numAngles;
347
    bool hasInitializedKernel;
348
    OpenCLContext& cl;
349
    const System& system;
350
    OpenCLArray* params;
351
352
};

353
354
355
356
357
/**
 * This kernel is invoked by CustomAngleForce to calculate the forces acting on the system and the energy of the system.
 */
class OpenCLCalcCustomAngleForceKernel : public CalcCustomAngleForceKernel {
public:
358
    OpenCLCalcCustomAngleForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomAngleForceKernel(name, platform),
359
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
360
361
362
363
364
365
366
367
368
369
    }
    ~OpenCLCalcCustomAngleForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomAngleForce this kernel will be used for
     */
    void initialize(const System& system, const CustomAngleForce& force);
    /**
370
     * Execute the kernel to calculate the forces and/or energy.
371
     *
372
373
374
375
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
376
     */
377
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
378
379
380
381
382
383
384
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomAngleForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomAngleForce& force);
385
386
387
388
private:
    int numAngles;
    bool hasInitializedKernel;
    OpenCLContext& cl;
389
    const System& system;
390
    OpenCLParameterSet* params;
391
    OpenCLArray* globals;
392
393
394
395
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

396
397
398
399
400
/**
 * This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class OpenCLCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
public:
401
    OpenCLCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcPeriodicTorsionForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
402
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
403
404
405
406
407
408
409
410
411
412
    }
    ~OpenCLCalcPeriodicTorsionForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the PeriodicTorsionForce this kernel will be used for
     */
    void initialize(const System& system, const PeriodicTorsionForce& force);
    /**
413
     * Execute the kernel to calculate the forces and/or energy.
414
     *
415
416
417
418
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
419
     */
420
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
421
422
423
424
425
426
427
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the PeriodicTorsionForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force);
428
429
private:
    int numTorsions;
430
    bool hasInitializedKernel;
431
    OpenCLContext& cl;
432
    const System& system;
433
    OpenCLArray* params;
434
435
};

436
437
438
439
440
/**
 * This kernel is invoked by RBTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class OpenCLCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
public:
441
    OpenCLCalcRBTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcRBTorsionForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
442
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
443
444
445
446
447
448
449
450
451
452
    }
    ~OpenCLCalcRBTorsionForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the RBTorsionForce this kernel will be used for
     */
    void initialize(const System& system, const RBTorsionForce& force);
    /**
453
     * Execute the kernel to calculate the forces and/or energy.
454
     *
455
456
457
458
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
459
     */
460
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
461
462
463
464
465
466
467
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the RBTorsionForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const RBTorsionForce& force);
468
469
private:
    int numTorsions;
470
    bool hasInitializedKernel;
471
    OpenCLContext& cl;
472
    const System& system;
473
    OpenCLArray* params;
474
475
};

476
477
478
479
480
/**
 * This kernel is invoked by CMAPTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class OpenCLCalcCMAPTorsionForceKernel : public CalcCMAPTorsionForceKernel {
public:
481
    OpenCLCalcCMAPTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCMAPTorsionForceKernel(name, platform),
482
            hasInitializedKernel(false), cl(cl), system(system), coefficients(NULL), mapPositions(NULL), torsionMaps(NULL) {
483
484
485
486
487
488
489
490
491
492
    }
    ~OpenCLCalcCMAPTorsionForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CMAPTorsionForce this kernel will be used for
     */
    void initialize(const System& system, const CMAPTorsionForce& force);
    /**
493
     * Execute the kernel to calculate the forces and/or energy.
494
     *
495
496
497
498
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
499
     */
500
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
501
502
503
504
private:
    int numTorsions;
    bool hasInitializedKernel;
    OpenCLContext& cl;
505
    const System& system;
506
507
508
    OpenCLArray* coefficients;
    OpenCLArray* mapPositions;
    OpenCLArray* torsionMaps;
509
510
};

511
512
513
514
515
/**
 * This kernel is invoked by CustomTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class OpenCLCalcCustomTorsionForceKernel : public CalcCustomTorsionForceKernel {
public:
516
    OpenCLCalcCustomTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomTorsionForceKernel(name, platform),
517
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
518
519
520
521
522
523
524
525
526
527
    }
    ~OpenCLCalcCustomTorsionForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomTorsionForce this kernel will be used for
     */
    void initialize(const System& system, const CustomTorsionForce& force);
    /**
528
     * Execute the kernel to calculate the forces and/or energy.
529
     *
530
531
532
533
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
534
     */
535
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
536
537
538
539
540
541
542
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomTorsionForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomTorsionForce& force);
543
544
545
546
private:
    int numTorsions;
    bool hasInitializedKernel;
    OpenCLContext& cl;
547
    const System& system;
548
    OpenCLParameterSet* params;
549
    OpenCLArray* globals;
550
551
552
553
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

554
555
556
557
558
/**
 * This kernel is invoked by NonbondedForce to calculate the forces acting on the system.
 */
class OpenCLCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
public:
559
    OpenCLCalcNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcNonbondedForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
560
            hasInitializedKernel(false), cl(cl), sigmaEpsilon(NULL), exceptionParams(NULL), cosSinSums(NULL), pmeGrid(NULL),
561
            pmeGrid2(NULL), pmeBsplineModuliX(NULL), pmeBsplineModuliY(NULL), pmeBsplineModuliZ(NULL), pmeBsplineTheta(NULL),
562
            pmeAtomRange(NULL), pmeAtomGridIndex(NULL), sort(NULL), fft(NULL), pmeio(NULL) {
563
564
565
566
567
568
569
570
571
572
    }
    ~OpenCLCalcNonbondedForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the NonbondedForce this kernel will be used for
     */
    void initialize(const System& system, const NonbondedForce& force);
    /**
573
     * Execute the kernel to calculate the forces and/or energy.
574
     *
575
576
577
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
578
579
     * @param includeDirect  true if direct space interactions should be included
     * @param includeReciprocal  true if reciprocal space interactions should be included
580
     * @return the potential energy due to the force
581
     */
582
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal);
583
584
585
586
587
588
589
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the NonbondedForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const NonbondedForce& force);
590
private:
591
592
593
594
595
596
597
598
599
    class SortTrait : public OpenCLSort::SortTrait {
        int getDataSize() const {return 8;}
        int getKeySize() const {return 4;}
        const char* getDataType() const {return "int2";}
        const char* getKeyType() const {return "int";}
        const char* getMinKey() const {return "INT_MIN";}
        const char* getMaxKey() const {return "INT_MAX";}
        const char* getMaxValue() const {return "(int2) (INT_MAX, INT_MAX)";}
        const char* getSortKey() const {return "value.y";}
600
    };
601
602
603
    class PmeIO;
    class PmePreComputation;
    class PmePostComputation;
604
605
    class SyncQueuePreComputation;
    class SyncQueuePostComputation;
606
    OpenCLContext& cl;
607
    bool hasInitializedKernel;
608
609
610
611
612
613
614
615
616
617
618
    OpenCLArray* sigmaEpsilon;
    OpenCLArray* exceptionParams;
    OpenCLArray* cosSinSums;
    OpenCLArray* pmeGrid;
    OpenCLArray* pmeGrid2;
    OpenCLArray* pmeBsplineModuliX;
    OpenCLArray* pmeBsplineModuliY;
    OpenCLArray* pmeBsplineModuliZ;
    OpenCLArray* pmeBsplineTheta;
    OpenCLArray* pmeAtomRange;
    OpenCLArray* pmeAtomGridIndex;
619
    OpenCLSort* sort;
620
621
    cl::CommandQueue pmeQueue;
    cl::Event pmeSyncEvent;
622
    OpenCLFFT3D* fft;
623
624
    Kernel cpuPme;
    PmeIO* pmeio;
625
626
    cl::Kernel ewaldSumsKernel;
    cl::Kernel ewaldForcesKernel;
627
628
    cl::Kernel pmeGridIndexKernel;
    cl::Kernel pmeAtomRangeKernel;
629
    cl::Kernel pmeZIndexKernel;
630
631
    cl::Kernel pmeUpdateBsplinesKernel;
    cl::Kernel pmeSpreadChargeKernel;
632
    cl::Kernel pmeFinishSpreadChargeKernel;
633
634
635
    cl::Kernel pmeConvolutionKernel;
    cl::Kernel pmeInterpolateForceKernel;
    std::map<std::string, std::string> pmeDefines;
636
    std::vector<std::pair<int, int> > exceptionAtoms;
637
    double ewaldSelfEnergy, dispersionCoefficient, alpha;
638
    bool hasCoulomb, hasLJ, usePmeQueue;
639
    static const int PmeOrder = 5;
640
641
};

642
643
644
645
646
/**
 * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
public:
647
    OpenCLCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomNonbondedForceKernel(name, platform),
648
            cl(cl), params(NULL), globals(NULL), interactionGroupData(NULL), forceCopy(NULL), system(system), hasInitializedKernel(false) {
649
650
651
652
653
654
655
656
657
658
    }
    ~OpenCLCalcCustomNonbondedForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomNonbondedForce this kernel will be used for
     */
    void initialize(const System& system, const CustomNonbondedForce& force);
    /**
659
     * Execute the kernel to calculate the forces and/or energy.
660
     *
661
662
663
664
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
665
     */
666
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
667
668
669
670
671
672
673
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomNonbondedForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
674
private:
675
    void initInteractionGroups(const CustomNonbondedForce& force, const std::string& interactionSource);
676
    OpenCLContext& cl;
677
    OpenCLParameterSet* params;
678
    OpenCLArray* globals;
679
680
681
    OpenCLArray* interactionGroupData;
    cl::Kernel interactionGroupKernel;
    std::vector<void*> interactionGroupArgs;
682
683
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
684
    std::vector<OpenCLArray*> tabulatedFunctions;
685
    double longRangeCoefficient;
686
687
    bool hasInitializedLongRangeCorrection, hasInitializedKernel;
    int numGroupThreadBlocks;
688
    CustomNonbondedForce* forceCopy;
689
    const System& system;
690
};
691
692
693
694
695
696

/**
 * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
 */
class OpenCLCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
public:
697
    OpenCLCalcGBSAOBCForceKernel(std::string name, const Platform& platform, OpenCLContext& cl) : CalcGBSAOBCForceKernel(name, platform), cl(cl),
698
699
            hasCreatedKernels(false), params(NULL), bornSum(NULL), longBornSum(NULL), bornRadii(NULL), bornForce(NULL),
            longBornForce(NULL), obcChain(NULL) {
700
701
702
703
704
705
706
707
708
709
    }
    ~OpenCLCalcGBSAOBCForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the GBSAOBCForce this kernel will be used for
     */
    void initialize(const System& system, const GBSAOBCForce& force);
    /**
710
     * Execute the kernel to calculate the forces and/or energy.
711
     *
712
713
714
715
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
716
     */
717
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
718
719
720
721
722
723
724
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the GBSAOBCForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
725
private:
726
    double prefactor, surfaceAreaFactor;
727
    bool hasCreatedKernels;
728
    int maxTiles;
729
    OpenCLContext& cl;
730
731
732
733
734
735
736
    OpenCLArray* params;
    OpenCLArray* bornSum;
    OpenCLArray* longBornSum;
    OpenCLArray* bornRadii;
    OpenCLArray* bornForce;
    OpenCLArray* longBornForce;
    OpenCLArray* obcChain;
737
738
    cl::Kernel computeBornSumKernel;
    cl::Kernel reduceBornSumKernel;
739
740
    cl::Kernel force1Kernel;
    cl::Kernel reduceBornForceKernel;
741
};
742

743
744
745
746
747
/**
 * This kernel is invoked by CustomGBForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomGBForceKernel : public CalcCustomGBForceKernel {
public:
748
    OpenCLCalcCustomGBForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomGBForceKernel(name, platform),
749
            hasInitializedKernels(false), cl(cl), params(NULL), computedValues(NULL), energyDerivs(NULL), energyDerivChain(NULL), longEnergyDerivs(NULL), globals(NULL),
750
            valueBuffers(NULL), longValueBuffers(NULL), system(system) {
751
752
753
754
755
756
757
758
759
760
    }
    ~OpenCLCalcCustomGBForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomGBForce this kernel will be used for
     */
    void initialize(const System& system, const CustomGBForce& force);
    /**
761
     * Execute the kernel to calculate the forces and/or energy.
762
     *
763
764
765
766
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
767
     */
768
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
769
770
771
772
773
774
775
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomGBForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomGBForce& force);
776
private:
Peter Eastman's avatar
Peter Eastman committed
777
    bool hasInitializedKernels, needParameterGradient;
778
    int maxTiles, numComputedValues;
779
780
781
    OpenCLContext& cl;
    OpenCLParameterSet* params;
    OpenCLParameterSet* computedValues;
782
    OpenCLParameterSet* energyDerivs;
783
    OpenCLParameterSet* energyDerivChain;
784
785
786
787
    OpenCLArray* longEnergyDerivs;
    OpenCLArray* globals;
    OpenCLArray* valueBuffers;
    OpenCLArray* longValueBuffers;
788
789
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
790
    std::vector<OpenCLArray*> tabulatedFunctions;
Peter Eastman's avatar
Peter Eastman committed
791
    std::vector<bool> pairValueUsesParam, pairEnergyUsesParam, pairEnergyUsesValue;
792
    const System& system;
Peter Eastman's avatar
Peter Eastman committed
793
    cl::Kernel pairValueKernel, perParticleValueKernel, pairEnergyKernel, perParticleEnergyKernel, gradientChainRuleKernel;
794
795
    std::string pairValueSrc, pairEnergySrc;
    std::map<std::string, std::string> pairValueDefines, pairEnergyDefines;
796
797
};

798
799
800
801
802
/**
 * This kernel is invoked by CustomExternalForce to calculate the forces acting on the system and the energy of the system.
 */
class OpenCLCalcCustomExternalForceKernel : public CalcCustomExternalForceKernel {
public:
803
    OpenCLCalcCustomExternalForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomExternalForceKernel(name, platform),
804
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
805
806
807
808
809
810
811
812
813
814
    }
    ~OpenCLCalcCustomExternalForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomExternalForce this kernel will be used for
     */
    void initialize(const System& system, const CustomExternalForce& force);
    /**
815
     * Execute the kernel to calculate the forces and/or energy.
816
     *
817
818
819
820
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
821
     */
822
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
823
824
825
826
827
828
829
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomExternalForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force);
830
831
832
833
private:
    int numParticles;
    bool hasInitializedKernel;
    OpenCLContext& cl;
834
    const System& system;
835
    OpenCLParameterSet* params;
836
    OpenCLArray* globals;
837
838
839
840
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

841
842
843
844
845
/**
 * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomHbondForceKernel : public CalcCustomHbondForceKernel {
public:
846
    OpenCLCalcCustomHbondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomHbondForceKernel(name, platform),
847
            hasInitializedKernel(false), cl(cl), donorParams(NULL), acceptorParams(NULL), donors(NULL), acceptors(NULL),
848
            donorBufferIndices(NULL), acceptorBufferIndices(NULL), globals(NULL), donorExclusions(NULL), acceptorExclusions(NULL), system(system) {
849
850
851
852
853
854
855
856
857
858
    }
    ~OpenCLCalcCustomHbondForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomHbondForce this kernel will be used for
     */
    void initialize(const System& system, const CustomHbondForce& force);
    /**
859
     * Execute the kernel to calculate the forces and/or energy.
860
     *
861
862
863
864
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
865
     */
866
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
867
868
869
870
871
872
873
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomHbondForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomHbondForce& force);
874
875
876
877
878
879
private:
    int numDonors, numAcceptors;
    bool hasInitializedKernel;
    OpenCLContext& cl;
    OpenCLParameterSet* donorParams;
    OpenCLParameterSet* acceptorParams;
880
881
882
883
884
885
886
    OpenCLArray* globals;
    OpenCLArray* donors;
    OpenCLArray* acceptors;
    OpenCLArray* donorBufferIndices;
    OpenCLArray* acceptorBufferIndices;
    OpenCLArray* donorExclusions;
    OpenCLArray* acceptorExclusions;
887
888
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
889
    std::vector<OpenCLArray*> tabulatedFunctions;
890
    const System& system;
891
    cl::Kernel donorKernel, acceptorKernel;
892
893
};

894
895
896
897
898
/**
 * This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomCompoundBondForceKernel : public CalcCustomCompoundBondForceKernel {
public:
899
    OpenCLCalcCustomCompoundBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomCompoundBondForceKernel(name, platform),
900
            cl(cl), params(NULL), globals(NULL), system(system) {
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
    }
    ~OpenCLCalcCustomCompoundBondForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomCompoundBondForce this kernel will be used for
     */
    void initialize(const System& system, const CustomCompoundBondForce& force);
    /**
     * Execute the kernel to calculate the forces and/or energy.
     *
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
     */
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
919
920
921
922
923
924
925
926
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomCompoundBondForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force);

927
928
929
930
private:
    int numBonds;
    OpenCLContext& cl;
    OpenCLParameterSet* params;
931
    OpenCLArray* globals;
932
933
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
934
    std::vector<OpenCLArray*> tabulatedFunctions;
935
    const System& system;
936
937
};

938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
/**
 * This kernel is invoked by CustomManyParticleForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomManyParticleForceKernel : public CalcCustomManyParticleForceKernel {
public:
    OpenCLCalcCustomManyParticleForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomManyParticleForceKernel(name, platform),
            hasInitializedKernel(false), cl(cl), params(NULL), globals(NULL), particleTypes(NULL), orderIndex(NULL), particleOrder(NULL), exclusions(NULL),
            exclusionStartIndex(NULL), blockCenter(NULL), blockBoundingBox(NULL), neighborPairs(NULL), numNeighborPairs(NULL), neighborStartIndex(NULL),
            numNeighborsForAtom(NULL), neighbors(NULL), system(system) {
    }
    ~OpenCLCalcCustomManyParticleForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomManyParticleForce this kernel will be used for
     */
    void initialize(const System& system, const CustomManyParticleForce& force);
    /**
     * Execute the kernel to calculate the forces and/or energy.
     *
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
     */
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomManyParticleForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force);

private:
    OpenCLContext& cl;
    bool hasInitializedKernel;
    NonbondedMethod nonbondedMethod;
    int maxNeighborPairs, forceWorkgroupSize, findNeighborsWorkgroupSize;
    OpenCLParameterSet* params;
    OpenCLArray* globals;
    OpenCLArray* particleTypes;
    OpenCLArray* orderIndex;
    OpenCLArray* particleOrder;
    OpenCLArray* exclusions;
    OpenCLArray* exclusionStartIndex;
    OpenCLArray* blockCenter;
    OpenCLArray* blockBoundingBox;
    OpenCLArray* neighborPairs;
    OpenCLArray* numNeighborPairs;
    OpenCLArray* neighborStartIndex;
    OpenCLArray* numNeighborsForAtom;
    OpenCLArray* neighbors;
    std::vector<std::string> globalParamNames;
    std::vector<float> globalParamValues;
    std::vector<OpenCLArray*> tabulatedFunctions;
    const System& system;
    cl::Kernel forceKernel, blockBoundsKernel, neighborsKernel, startIndicesKernel, copyPairsKernel;
};

999
1000
1001
1002
1003
/**
 * This kernel is invoked by VerletIntegrator to take one time step.
 */
class OpenCLIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
public:
1004
    OpenCLIntegrateVerletStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateVerletStepKernel(name, platform), cl(cl),
1005
            hasInitializedKernels(false) {
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
    }
    ~OpenCLIntegrateVerletStepKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param integrator the VerletIntegrator this kernel will be used for
     */
    void initialize(const System& system, const VerletIntegrator& integrator);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
     * @param integrator the VerletIntegrator this kernel is being used for
     */
    void execute(ContextImpl& context, const VerletIntegrator& integrator);
1022
1023
1024
1025
1026
1027
1028
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the VerletIntegrator this kernel is being used for
     */
    double computeKineticEnergy(ContextImpl& context, const VerletIntegrator& integrator);
1029
private:
1030
    OpenCLContext& cl;
1031
    double prevStepSize;
1032
    bool hasInitializedKernels;
1033
1034
1035
1036
1037
1038
1039
1040
1041
    cl::Kernel kernel1, kernel2;
};

/**
 * This kernel is invoked by LangevinIntegrator to take one time step.
 */
class OpenCLIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
public:
    OpenCLIntegrateLangevinStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateLangevinStepKernel(name, platform), cl(cl),
1042
            hasInitializedKernels(false), params(NULL) {
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
    }
    ~OpenCLIntegrateLangevinStepKernel();
    /**
     * Initialize the kernel, setting up the particle masses.
     *
     * @param system     the System this kernel will be applied to
     * @param integrator the LangevinIntegrator this kernel will be used for
     */
    void initialize(const System& system, const LangevinIntegrator& integrator);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
     * @param integrator the LangevinIntegrator this kernel is being used for
     */
    void execute(ContextImpl& context, const LangevinIntegrator& integrator);
1059
1060
1061
1062
1063
1064
1065
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the LangevinIntegrator this kernel is being used for
     */
    double computeKineticEnergy(ContextImpl& context, const LangevinIntegrator& integrator);
1066
1067
1068
private:
    OpenCLContext& cl;
    double prevTemp, prevFriction, prevStepSize;
1069
    bool hasInitializedKernels;
1070
    OpenCLArray* params;
1071
    cl::Kernel kernel1, kernel2;
1072
1073
};

1074
1075
1076
1077
1078
/**
 * This kernel is invoked by BrownianIntegrator to take one time step.
 */
class OpenCLIntegrateBrownianStepKernel : public IntegrateBrownianStepKernel {
public:
1079
1080
    OpenCLIntegrateBrownianStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateBrownianStepKernel(name, platform), cl(cl),
            hasInitializedKernels(false), prevTemp(-1), prevFriction(-1), prevStepSize(-1) {
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
    }
    ~OpenCLIntegrateBrownianStepKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param integrator the BrownianIntegrator this kernel will be used for
     */
    void initialize(const System& system, const BrownianIntegrator& integrator);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
     * @param integrator the BrownianIntegrator this kernel is being used for
     */
    void execute(ContextImpl& context, const BrownianIntegrator& integrator);
1097
1098
1099
1100
1101
1102
1103
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the BrownianIntegrator this kernel is being used for
     */
    double computeKineticEnergy(ContextImpl& context, const BrownianIntegrator& integrator);
1104
1105
1106
1107
1108
1109
private:
    OpenCLContext& cl;
    double prevTemp, prevFriction, prevStepSize;
    bool hasInitializedKernels;
    cl::Kernel kernel1, kernel2;
};
1110
1111
1112
1113
1114
1115
1116

/**
 * This kernel is invoked by VariableVerletIntegrator to take one time step.
 */
class OpenCLIntegrateVariableVerletStepKernel : public IntegrateVariableVerletStepKernel {
public:
    OpenCLIntegrateVariableVerletStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateVariableVerletStepKernel(name, platform), cl(cl),
1117
            hasInitializedKernels(false) {
1118
1119
1120
1121
1122
1123
    }
    ~OpenCLIntegrateVariableVerletStepKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
1124
     * @param integrator the VariableVerletIntegrator this kernel will be used for
1125
1126
1127
1128
1129
1130
     */
    void initialize(const System& system, const VariableVerletIntegrator& integrator);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
1131
     * @param integrator the VariableVerletIntegrator this kernel is being used for
1132
     * @param maxTime    the maximum time beyond which the simulation should not be advanced
1133
     * @return the size of the step that was taken
1134
     */
1135
    double execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
1136
1137
1138
1139
1140
1141
1142
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the VariableVerletIntegrator this kernel is being used for
     */
    double computeKineticEnergy(ContextImpl& context, const VariableVerletIntegrator& integrator);
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
private:
    OpenCLContext& cl;
    bool hasInitializedKernels;
    int blockSize;
    cl::Kernel kernel1, kernel2, selectSizeKernel;
};

/**
 * This kernel is invoked by VariableLangevinIntegrator to take one time step.
 */
class OpenCLIntegrateVariableLangevinStepKernel : public IntegrateVariableLangevinStepKernel {
public:
    OpenCLIntegrateVariableLangevinStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateVariableLangevinStepKernel(name, platform), cl(cl),
1156
            hasInitializedKernels(false), params(NULL) {
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
    }
    ~OpenCLIntegrateVariableLangevinStepKernel();
    /**
     * Initialize the kernel, setting up the particle masses.
     *
     * @param system     the System this kernel will be applied to
     * @param integrator the VariableLangevinIntegrator this kernel will be used for
     */
    void initialize(const System& system, const VariableLangevinIntegrator& integrator);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
     * @param integrator the VariableLangevinIntegrator this kernel is being used for
     * @param maxTime    the maximum time beyond which the simulation should not be advanced
1172
     * @return the size of the step that was taken
1173
     */
1174
    double execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
1175
1176
1177
1178
1179
1180
1181
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the VariableLangevinIntegrator this kernel is being used for
     */
    double computeKineticEnergy(ContextImpl& context, const VariableLangevinIntegrator& integrator);
1182
1183
1184
1185
private:
    OpenCLContext& cl;
    bool hasInitializedKernels;
    int blockSize;
1186
    OpenCLArray* params;
1187
    cl::Kernel kernel1, kernel2, selectSizeKernel;
1188
1189
    double prevTemp, prevFriction, prevErrorTol;
};
1190

1191
1192
1193
1194
1195
1196
/**
 * This kernel is invoked by CustomIntegrator to take one time step.
 */
class OpenCLIntegrateCustomStepKernel : public IntegrateCustomStepKernel {
public:
    OpenCLIntegrateCustomStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateCustomStepKernel(name, platform), cl(cl),
1197
1198
            hasInitializedKernels(false), localValuesAreCurrent(false), globalValues(NULL), contextParameterValues(NULL), sumBuffer(NULL), potentialEnergy(NULL),
            kineticEnergy(NULL), uniformRandoms(NULL), randomSeed(NULL), perDofValues(NULL) {
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
    }
    ~OpenCLIntegrateCustomStepKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param integrator the CustomIntegrator this kernel will be used for
     */
    void initialize(const System& system, const CustomIntegrator& integrator);
    /**
     * Execute the kernel.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the CustomIntegrator this kernel is being used for
     * @param forcesAreValid if the context has been modified since the last time step, this will be
     *                       false to show that cached forces are invalid and must be recalculated.
     *                       On exit, this should specify whether the cached forces are valid at the
     *                       end of the step.
     */
    void execute(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid);
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the CustomIntegrator this kernel is being used for
     * @param forcesAreValid if the context has been modified since the last time step, this will be
     *                       false to show that cached forces are invalid and must be recalculated.
     *                       On exit, this should specify whether the cached forces are valid at the
     *                       end of the step.
     */
    double computeKineticEnergy(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid);
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
    /**
     * Get the values of all global variables.
     *
     * @param context   the context in which to execute this kernel
     * @param values    on exit, this contains the values
     */
    void getGlobalVariables(ContextImpl& context, std::vector<double>& values) const;
    /**
     * Set the values of all global variables.
     *
     * @param context   the context in which to execute this kernel
     * @param values    a vector containing the values
     */
    void setGlobalVariables(ContextImpl& context, const std::vector<double>& values);
    /**
     * Get the values of a per-DOF variable.
     *
     * @param context   the context in which to execute this kernel
     * @param variable  the index of the variable to get
     * @param values    on exit, this contains the values
     */
    void getPerDofVariable(ContextImpl& context, int variable, std::vector<Vec3>& values) const;
    /**
     * Set the values of a per-DOF variable.
     *
     * @param context   the context in which to execute this kernel
     * @param variable  the index of the variable to get
     * @param values    a vector containing the values
     */
    void setPerDofVariable(ContextImpl& context, int variable, const std::vector<Vec3>& values);
private:
1261
    class ReorderListener;
1262
1263
    std::string createGlobalComputation(const std::string& variable, const Lepton::ParsedExpression& expr, CustomIntegrator& integrator, const std::string& energyName);
    std::string createPerDofComputation(const std::string& variable, const Lepton::ParsedExpression& expr, int component, CustomIntegrator& integrator, const std::string& forceName, const std::string& energyName);
1264
    void prepareForComputation(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid);
1265
    void recordChangedParameters(ContextImpl& context);
1266
    OpenCLContext& cl;
1267
    double prevStepSize;
1268
    int numGlobalVariables;
1269
    bool hasInitializedKernels, deviceValuesAreCurrent, modifiesParameters, keNeedsForce;
1270
    mutable bool localValuesAreCurrent;
1271
1272
1273
    OpenCLArray* globalValues;
    OpenCLArray* contextParameterValues;
    OpenCLArray* sumBuffer;
1274
1275
    OpenCLArray* potentialEnergy;
    OpenCLArray* kineticEnergy;
1276
1277
    OpenCLArray* uniformRandoms;
    OpenCLArray* randomSeed;
1278
1279
    std::map<int, OpenCLArray*> savedForces;
    std::set<int> validSavedForces;
1280
    OpenCLParameterSet* perDofValues;
1281
1282
1283
1284
    mutable std::vector<std::vector<cl_float> > localPerDofValuesFloat;
    mutable std::vector<std::vector<cl_double> > localPerDofValuesDouble;
    std::vector<float> contextValuesFloat;
    std::vector<double> contextValuesDouble;
1285
    std::vector<float> contextValues;
1286
    std::vector<std::vector<cl::Kernel> > kernels;
1287
    cl::Kernel sumPotentialEnergyKernel, randomKernel, kineticEnergyKernel, sumKineticEnergyKernel;
1288
1289
1290
1291
    std::vector<CustomIntegrator::ComputationType> stepType;
    std::vector<bool> needsForces;
    std::vector<bool> needsEnergy;
    std::vector<bool> invalidatesForces;
1292
    std::vector<bool> merged;
1293
    std::vector<int> forceGroup;
1294
1295
    std::vector<int> requiredGaussian;
    std::vector<int> requiredUniform;
1296
    std::vector<std::string> parameterNames;
1297
1298
};

1299
1300
1301
1302
1303
1304
/**
 * This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
 */
class OpenCLApplyAndersenThermostatKernel : public ApplyAndersenThermostatKernel {
public:
    OpenCLApplyAndersenThermostatKernel(std::string name, const Platform& platform, OpenCLContext& cl) : ApplyAndersenThermostatKernel(name, platform), cl(cl),
1305
            hasInitializedKernels(false), atomGroups(NULL) {
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
    }
    ~OpenCLApplyAndersenThermostatKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param thermostat the AndersenThermostat this kernel will be used for
     */
    void initialize(const System& system, const AndersenThermostat& thermostat);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
     */
    void execute(ContextImpl& context);
private:
    OpenCLContext& cl;
    bool hasInitializedKernels;
    int randomSeed;
1325
    OpenCLArray* atomGroups;
1326
    cl::Kernel kernel;
1327
1328
1329
1330
1331
1332
1333
1334
};

/**
 * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
 */
class OpenCLApplyMonteCarloBarostatKernel : public ApplyMonteCarloBarostatKernel {
public:
    OpenCLApplyMonteCarloBarostatKernel(std::string name, const Platform& platform, OpenCLContext& cl) : ApplyMonteCarloBarostatKernel(name, platform), cl(cl),
1335
            hasInitializedKernels(false), savedPositions(NULL), moleculeAtoms(NULL), moleculeStartIndex(NULL) {
1336
1337
1338
1339
1340
1341
1342
1343
    }
    ~OpenCLApplyMonteCarloBarostatKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param barostat   the MonteCarloBarostat this kernel will be used for
     */
1344
    void initialize(const System& system, const Force& barostat);
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
    /**
     * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
     * This version scales the x, y, and z positions independently.
     * This is called BEFORE the periodic box size is modified.  It should begin by translating each particle
     * or cluster into the first periodic box, so that coordinates will still be correct after the box size
     * is changed.
     *
     * @param context    the context in which to execute this kernel
     * @param scaleX     the scale factor by which to multiply particle x-coordinate
     * @param scaleY     the scale factor by which to multiply particle y-coordinate
     * @param scaleZ     the scale factor by which to multiply particle z-coordinate
     */
1357
    void scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ);
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
    /**
     * Reject the most recent Monte Carlo step, restoring the particle positions to where they were before
     * scaleCoordinates() was last called.
     *
     * @param context    the context in which to execute this kernel
     */
    void restoreCoordinates(ContextImpl& context);
private:
    OpenCLContext& cl;
    bool hasInitializedKernels;
    int numMolecules;
1369
1370
1371
    OpenCLArray* savedPositions;
    OpenCLArray* moleculeAtoms;
    OpenCLArray* moleculeStartIndex;
1372
    cl::Kernel kernel;
1373
    std::vector<int> lastAtomOrder;
1374
};
1375

1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
/**
 * This kernel is invoked to remove center of mass motion from the system.
 */
class OpenCLRemoveCMMotionKernel : public RemoveCMMotionKernel {
public:
    OpenCLRemoveCMMotionKernel(std::string name, const Platform& platform, OpenCLContext& cl) : RemoveCMMotionKernel(name, platform), cl(cl), cmMomentum(NULL) {
    }
    ~OpenCLRemoveCMMotionKernel();
    /**
     * Initialize the kernel, setting up the particle masses.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CMMotionRemover this kernel will be used for
     */
    void initialize(const System& system, const CMMotionRemover& force);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
     */
    void execute(ContextImpl& context);
private:
    OpenCLContext& cl;
    int frequency;
1400
    OpenCLArray* cmMomentum;
1401
1402
    cl::Kernel kernel1, kernel2;
};
1403
1404
1405
1406

} // namespace OpenMM

#endif /*OPENMM_OPENCLKERNELS_H_*/