OpenCLKernels.h 51.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-2010 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
     * @return the potential energy of the system.  This value is added to all values returned by ForceImpls'
75
     * calcForcesAndEnergy() methods.  That is, each force kernel may <i>either</i> return its contribution to the
76
77
     * energy directly, <i>or</i> add it to an internal buffer so that it will be included here.
     */
78
    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
79
private:
80
   OpenCLContext& cl;
81
82
83
};

/**
84
85
 * This kernel provides methods for setting and retrieving various state data: time, positions,
 * velocities, and forces.
86
 */
87
class OpenCLUpdateStateDataKernel : public UpdateStateDataKernel {
88
public:
89
    OpenCLUpdateStateDataKernel(std::string name, const Platform& platform, OpenCLContext& cl) : UpdateStateDataKernel(name, platform), cl(cl) {
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
    }
    /**
     * 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);
109
110
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
    /**
     * 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);
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
    /**
     * 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
155
156
157
158
159
160
161
162
163
164
165
166
    /**
     * 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);
167
private:
168
    OpenCLContext& cl;
169
};
170

171
172
173
174
175
/**
 * This kernel modifies the positions of particles to enforce distance constraints.
 */
class OpenCLApplyConstraintsKernel : public ApplyConstraintsKernel {
public:
176
177
    OpenCLApplyConstraintsKernel(std::string name, const Platform& platform, OpenCLContext& cl) : ApplyConstraintsKernel(name, platform),
            cl(cl), hasInitializedKernel(false) {
178
179
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);
private:
    OpenCLContext& cl;
194
195
    bool hasInitializedKernel;
    cl::Kernel applyDeltasKernel;
196
197
};

198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
/**
 * 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;
};

221
222
223
224
225
/**
 * 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:
226
    OpenCLCalcHarmonicBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcHarmonicBondForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
227
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
228
229
230
231
232
233
234
235
236
237
    }
    ~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);
    /**
238
     * Execute the kernel to calculate the forces and/or energy.
239
     *
240
241
242
243
     * @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
244
     */
245
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
246
247
248
249
250
251
252
    /**
     * 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);
253
254
private:
    int numBonds;
255
    bool hasInitializedKernel;
256
    OpenCLContext& cl;
257
258
259
260
    System& system;
    OpenCLArray<mm_float2>* params;
};

261
262
263
264
265
266
/**
 * 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:
    OpenCLCalcCustomBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomBondForceKernel(name, platform),
267
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
268
269
270
271
272
273
274
275
276
277
    }
    ~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);
    /**
278
     * Execute the kernel to calculate the forces and/or energy.
279
     *
280
281
282
283
     * @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
284
     */
285
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
286
287
288
289
290
291
292
    /**
     * 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);
293
294
295
296
297
private:
    int numBonds;
    bool hasInitializedKernel;
    OpenCLContext& cl;
    System& system;
298
    OpenCLParameterSet* params;
299
300
301
302
303
    OpenCLArray<cl_float>* globals;
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

304
305
306
307
308
/**
 * 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:
309
    OpenCLCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcHarmonicAngleForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
310
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
311
312
313
314
315
316
317
318
319
320
    }
    ~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);
    /**
321
     * Execute the kernel to calculate the forces and/or energy.
322
     *
323
324
325
326
     * @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
327
     */
328
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
329
330
331
332
333
334
335
    /**
     * 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);
336
337
private:
    int numAngles;
338
    bool hasInitializedKernel;
339
    OpenCLContext& cl;
340
341
342
343
    System& system;
    OpenCLArray<mm_float2>* params;
};

344
345
346
347
348
349
/**
 * 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:
    OpenCLCalcCustomAngleForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomAngleForceKernel(name, platform),
350
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
351
352
353
354
355
356
357
358
359
360
    }
    ~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);
    /**
361
     * Execute the kernel to calculate the forces and/or energy.
362
     *
363
364
365
366
     * @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
367
     */
368
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
369
370
371
372
373
374
375
    /**
     * 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);
376
377
378
379
380
381
382
383
384
385
386
private:
    int numAngles;
    bool hasInitializedKernel;
    OpenCLContext& cl;
    System& system;
    OpenCLParameterSet* params;
    OpenCLArray<cl_float>* globals;
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

387
388
389
390
391
/**
 * 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:
392
    OpenCLCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcPeriodicTorsionForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
393
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
394
395
396
397
398
399
400
401
402
403
    }
    ~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);
    /**
404
     * Execute the kernel to calculate the forces and/or energy.
405
     *
406
407
408
409
     * @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
410
     */
411
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
412
413
414
415
416
417
418
    /**
     * 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);
419
420
private:
    int numTorsions;
421
    bool hasInitializedKernel;
422
    OpenCLContext& cl;
423
424
425
426
    System& system;
    OpenCLArray<mm_float4>* params;
};

427
428
429
430
431
/**
 * 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:
432
    OpenCLCalcRBTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcRBTorsionForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
433
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
434
435
436
437
438
439
440
441
442
443
    }
    ~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);
    /**
444
     * Execute the kernel to calculate the forces and/or energy.
445
     *
446
447
448
449
     * @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
450
     */
451
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
452
453
454
455
456
457
458
    /**
     * 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);
459
460
private:
    int numTorsions;
461
    bool hasInitializedKernel;
462
    OpenCLContext& cl;
463
464
465
466
    System& system;
    OpenCLArray<mm_float8>* params;
};

467
468
469
470
471
472
/**
 * 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:
    OpenCLCalcCMAPTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCMAPTorsionForceKernel(name, platform),
473
            hasInitializedKernel(false), cl(cl), system(system), coefficients(NULL), mapPositions(NULL), torsionMaps(NULL) {
474
475
476
477
478
479
480
481
482
483
    }
    ~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);
    /**
484
     * Execute the kernel to calculate the forces and/or energy.
485
     *
486
487
488
489
     * @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
490
     */
491
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
492
493
494
495
496
497
498
499
500
501
private:
    int numTorsions;
    bool hasInitializedKernel;
    OpenCLContext& cl;
    System& system;
    OpenCLArray<mm_float4>* coefficients;
    OpenCLArray<mm_int2>* mapPositions;
    OpenCLArray<cl_int>* torsionMaps;
};

502
503
504
505
506
507
/**
 * 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:
    OpenCLCalcCustomTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomTorsionForceKernel(name, platform),
508
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
509
510
511
512
513
514
515
516
517
518
    }
    ~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);
    /**
519
     * Execute the kernel to calculate the forces and/or energy.
520
     *
521
522
523
524
     * @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
525
     */
526
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
527
528
529
530
531
532
533
    /**
     * 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);
534
535
536
537
538
539
540
541
542
543
544
private:
    int numTorsions;
    bool hasInitializedKernel;
    OpenCLContext& cl;
    System& system;
    OpenCLParameterSet* params;
    OpenCLArray<cl_float>* globals;
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

545
546
547
548
549
/**
 * This kernel is invoked by NonbondedForce to calculate the forces acting on the system.
 */
class OpenCLCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
public:
550
    OpenCLCalcNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcNonbondedForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
551
            hasInitializedKernel(false), cl(cl), sigmaEpsilon(NULL), exceptionParams(NULL), cosSinSums(NULL), pmeGrid(NULL),
552
553
            pmeGrid2(NULL), pmeBsplineModuliX(NULL), pmeBsplineModuliY(NULL), pmeBsplineModuliZ(NULL), pmeBsplineTheta(NULL), pmeBsplineDTheta(NULL),
            pmeAtomRange(NULL), pmeAtomGridIndex(NULL), sort(NULL), fft(NULL) {
554
555
556
557
558
559
560
561
562
563
    }
    ~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);
    /**
564
     * Execute the kernel to calculate the forces and/or energy.
565
     *
566
567
568
     * @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
569
570
     * @param includeDirect  true if direct space interactions should be included
     * @param includeReciprocal  true if reciprocal space interactions should be included
571
     * @return the potential energy due to the force
572
     */
573
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal);
574
575
576
577
578
579
580
    /**
     * 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);
581
private:
582
583
584
585
586
587
588
589
590
591
    struct SortTrait {
        typedef mm_int2 DataType;
        typedef cl_int KeyType;
        static const char* clDataType() {return "int2";}
        static const char* clKeyType() {return "int";}
        static const char* clMinKey() {return "INT_MIN";}
        static const char* clMaxKey() {return "INT_MAX";}
        static const char* clMaxValue() {return "(int2) (INT_MAX, INT_MAX)";}
        static const char* clSortKey() {return "value.y";}
    };
592
    OpenCLContext& cl;
593
    bool hasInitializedKernel;
594
    OpenCLArray<mm_float2>* sigmaEpsilon;
595
    OpenCLArray<mm_float4>* exceptionParams;
596
    OpenCLArray<mm_float2>* cosSinSums;
597
    OpenCLArray<mm_float2>* pmeGrid;
Peter Eastman's avatar
Peter Eastman committed
598
    OpenCLArray<mm_float2>* pmeGrid2;
599
600
601
602
    OpenCLArray<cl_float>* pmeBsplineModuliX;
    OpenCLArray<cl_float>* pmeBsplineModuliY;
    OpenCLArray<cl_float>* pmeBsplineModuliZ;
    OpenCLArray<mm_float4>* pmeBsplineTheta;
603
    OpenCLArray<mm_float4>* pmeBsplineDTheta;
604
    OpenCLArray<cl_int>* pmeAtomRange;
605
    OpenCLArray<mm_int2>* pmeAtomGridIndex;
606
    OpenCLSort<SortTrait>* sort;
607
    OpenCLFFT3D* fft;
608
609
    cl::Kernel ewaldSumsKernel;
    cl::Kernel ewaldForcesKernel;
610
611
    cl::Kernel pmeGridIndexKernel;
    cl::Kernel pmeAtomRangeKernel;
612
    cl::Kernel pmeZIndexKernel;
613
614
    cl::Kernel pmeUpdateBsplinesKernel;
    cl::Kernel pmeSpreadChargeKernel;
615
    cl::Kernel pmeFinishSpreadChargeKernel;
616
617
618
    cl::Kernel pmeConvolutionKernel;
    cl::Kernel pmeInterpolateForceKernel;
    std::map<std::string, std::string> pmeDefines;
619
    double ewaldSelfEnergy, dispersionCoefficient, alpha;
620
    int interpolateForceThreads;
621
    bool hasCoulomb, hasLJ;
622
    static const int PmeOrder = 5;
623
624
};

625
626
627
628
629
630
/**
 * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
public:
    OpenCLCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomNonbondedForceKernel(name, platform),
631
            hasInitializedKernel(false), cl(cl), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), system(system) {
632
633
634
635
636
637
638
639
640
641
    }
    ~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);
    /**
642
     * Execute the kernel to calculate the forces and/or energy.
643
     *
644
645
646
647
     * @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
648
     */
649
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
650
private:
651
    bool hasInitializedKernel;
652
    OpenCLContext& cl;
653
    OpenCLParameterSet* params;
654
    OpenCLArray<cl_float>* globals;
655
    OpenCLArray<mm_float4>* tabulatedFunctionParams;
656
657
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
658
    std::vector<OpenCLArray<mm_float4>*> tabulatedFunctions;
659
660
    System& system;
};
661
662
663
664
665
666

/**
 * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
 */
class OpenCLCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
public:
667
    OpenCLCalcGBSAOBCForceKernel(std::string name, const Platform& platform, OpenCLContext& cl) : CalcGBSAOBCForceKernel(name, platform), cl(cl),
668
669
            hasCreatedKernels(false), params(NULL), bornSum(NULL), longBornSum(NULL), bornRadii(NULL), bornForce(NULL),
            longBornForce(NULL), obcChain(NULL) {
670
671
672
673
674
675
676
677
678
679
    }
    ~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);
    /**
680
     * Execute the kernel to calculate the forces and/or energy.
681
     *
682
683
684
685
     * @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
686
     */
687
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
688
689
690
691
692
693
694
    /**
     * 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);
695
696
private:
    double prefactor;
697
    bool hasCreatedKernels;
698
    int maxTiles;
699
700
701
    OpenCLContext& cl;
    OpenCLArray<mm_float2>* params;
    OpenCLArray<cl_float>* bornSum;
702
    OpenCLArray<cl_long>* longBornSum;
703
704
    OpenCLArray<cl_float>* bornRadii;
    OpenCLArray<cl_float>* bornForce;
705
    OpenCLArray<cl_long>* longBornForce;
706
707
708
    OpenCLArray<cl_float>* obcChain;
    cl::Kernel computeBornSumKernel;
    cl::Kernel reduceBornSumKernel;
709
710
    cl::Kernel force1Kernel;
    cl::Kernel reduceBornForceKernel;
711
};
712

713
714
715
716
717
718
/**
 * This kernel is invoked by CustomGBForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomGBForceKernel : public CalcCustomGBForceKernel {
public:
    OpenCLCalcCustomGBForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomGBForceKernel(name, platform),
719
720
            hasInitializedKernels(false), cl(cl), params(NULL), computedValues(NULL), energyDerivs(NULL), longEnergyDerivs(NULL), globals(NULL),
            valueBuffers(NULL), longValueBuffers(NULL), tabulatedFunctionParams(NULL), system(system) {
721
722
723
724
725
726
727
728
729
730
    }
    ~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);
    /**
731
     * Execute the kernel to calculate the forces and/or energy.
732
     *
733
734
735
736
     * @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
737
     */
738
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
739
private:
Peter Eastman's avatar
Peter Eastman committed
740
    bool hasInitializedKernels, needParameterGradient;
741
    int maxTiles, numComputedValues;
742
743
744
    OpenCLContext& cl;
    OpenCLParameterSet* params;
    OpenCLParameterSet* computedValues;
745
    OpenCLParameterSet* energyDerivs;
746
    OpenCLArray<cl_long>* longEnergyDerivs;
747
748
    OpenCLArray<cl_float>* globals;
    OpenCLArray<cl_float>* valueBuffers;
749
    OpenCLArray<cl_long>* longValueBuffers;
750
751
752
753
    OpenCLArray<mm_float4>* tabulatedFunctionParams;
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
    std::vector<OpenCLArray<mm_float4>*> tabulatedFunctions;
Peter Eastman's avatar
Peter Eastman committed
754
    std::vector<bool> pairValueUsesParam, pairEnergyUsesParam, pairEnergyUsesValue;
755
    System& system;
Peter Eastman's avatar
Peter Eastman committed
756
    cl::Kernel pairValueKernel, perParticleValueKernel, pairEnergyKernel, perParticleEnergyKernel, gradientChainRuleKernel;
757
758
};

759
760
761
762
763
764
/**
 * 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:
    OpenCLCalcCustomExternalForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomExternalForceKernel(name, platform),
765
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
766
767
768
769
770
771
772
773
774
775
    }
    ~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);
    /**
776
     * Execute the kernel to calculate the forces and/or energy.
777
     *
778
779
780
781
     * @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
782
     */
783
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
784
785
786
787
788
private:
    int numParticles;
    bool hasInitializedKernel;
    OpenCLContext& cl;
    System& system;
789
    OpenCLParameterSet* params;
790
791
792
793
794
    OpenCLArray<cl_float>* globals;
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

795
796
797
798
799
800
801
/**
 * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomHbondForceKernel : public CalcCustomHbondForceKernel {
public:
    OpenCLCalcCustomHbondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomHbondForceKernel(name, platform),
            hasInitializedKernel(false), cl(cl), donorParams(NULL), acceptorParams(NULL), donors(NULL), acceptors(NULL),
802
803
            donorBufferIndices(NULL), acceptorBufferIndices(NULL), globals(NULL), donorExclusions(NULL), acceptorExclusions(NULL),
            tabulatedFunctionParams(NULL), system(system) {
804
805
806
807
808
809
810
811
812
813
    }
    ~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);
    /**
814
     * Execute the kernel to calculate the forces and/or energy.
815
     *
816
817
818
819
     * @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
820
     */
821
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
822
823
824
825
826
827
828
829
830
831
832
private:
    int numDonors, numAcceptors;
    bool hasInitializedKernel;
    OpenCLContext& cl;
    OpenCLParameterSet* donorParams;
    OpenCLParameterSet* acceptorParams;
    OpenCLArray<cl_float>* globals;
    OpenCLArray<mm_int4>* donors;
    OpenCLArray<mm_int4>* acceptors;
    OpenCLArray<mm_int4>* donorBufferIndices;
    OpenCLArray<mm_int4>* acceptorBufferIndices;
833
834
    OpenCLArray<mm_int4>* donorExclusions;
    OpenCLArray<mm_int4>* acceptorExclusions;
835
836
837
838
839
    OpenCLArray<mm_float4>* tabulatedFunctionParams;
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
    std::vector<OpenCLArray<mm_float4>*> tabulatedFunctions;
    System& system;
840
    cl::Kernel donorKernel, acceptorKernel;
841
842
};

843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
/**
 * This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomCompoundBondForceKernel : public CalcCustomCompoundBondForceKernel {
public:
    OpenCLCalcCustomCompoundBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomCompoundBondForceKernel(name, platform),
            cl(cl), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), system(system) {
    }
    ~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);
private:
    int numBonds;
    OpenCLContext& cl;
    OpenCLParameterSet* params;
    OpenCLArray<cl_float>* globals;
    OpenCLArray<mm_float4>* tabulatedFunctionParams;
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
    std::vector<OpenCLArray<mm_float4>*> tabulatedFunctions;
    System& system;
};

880
881
882
883
884
/**
 * This kernel is invoked by VerletIntegrator to take one time step.
 */
class OpenCLIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
public:
885
    OpenCLIntegrateVerletStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateVerletStepKernel(name, platform), cl(cl),
886
            hasInitializedKernels(false) {
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
    }
    ~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);
private:
904
    OpenCLContext& cl;
905
    double prevStepSize;
906
    bool hasInitializedKernels;
907
908
909
910
911
912
913
914
915
    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),
916
            hasInitializedKernels(false), params(NULL) {
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
    }
    ~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);
private:
    OpenCLContext& cl;
    double prevTemp, prevFriction, prevStepSize;
936
    bool hasInitializedKernels;
937
    OpenCLArray<cl_float>* params;
938
    cl::Kernel kernel1, kernel2;
939
940
};

941
942
943
944
945
/**
 * This kernel is invoked by BrownianIntegrator to take one time step.
 */
class OpenCLIntegrateBrownianStepKernel : public IntegrateBrownianStepKernel {
public:
946
947
    OpenCLIntegrateBrownianStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateBrownianStepKernel(name, platform), cl(cl),
            hasInitializedKernels(false), prevTemp(-1), prevFriction(-1), prevStepSize(-1) {
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
    }
    ~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);
private:
    OpenCLContext& cl;
    double prevTemp, prevFriction, prevStepSize;
    bool hasInitializedKernels;
    cl::Kernel kernel1, kernel2;
};
970
971
972
973
974
975
976

/**
 * 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),
977
            hasInitializedKernels(false) {
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
    }
    ~OpenCLIntegrateVariableVerletStepKernel();
    /**
     * 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 VariableVerletIntegrator& integrator);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
     * @param integrator the VerletIntegrator this kernel is being used for
     * @param maxTime    the maximum time beyond which the simulation should not be advanced
993
     * @return the size of the step that was taken
994
     */
995
    double execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
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),
1009
            hasInitializedKernels(false), params(NULL) {
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
    }
    ~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
1025
     * @return the size of the step that was taken
1026
     */
1027
    double execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
1028
1029
1030
1031
1032
private:
    OpenCLContext& cl;
    bool hasInitializedKernels;
    int blockSize;
    OpenCLArray<cl_float>* params;
1033
    cl::Kernel kernel1, kernel2, selectSizeKernel;
1034
1035
    double prevTemp, prevFriction, prevErrorTol;
};
1036

1037
1038
1039
1040
1041
1042
/**
 * 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),
1043
1044
            hasInitializedKernels(false), localValuesAreCurrent(false), globalValues(NULL), contextParameterValues(NULL), sumBuffer(NULL), energy(NULL),
            uniformRandoms(NULL), randomSeed(NULL), perDofValues(NULL) {
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
    }
    ~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);
    /**
     * 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:
1096
    class ReorderListener;
1097
1098
    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);
1099
    void recordChangedParameters(ContextImpl& context);
1100
    OpenCLContext& cl;
1101
    double prevStepSize;
1102
    int numGlobalVariables;
1103
    bool hasInitializedKernels, deviceValuesAreCurrent, modifiesParameters;
1104
1105
    mutable bool localValuesAreCurrent;
    OpenCLArray<cl_float>* globalValues;
1106
1107
    OpenCLArray<cl_float>* contextParameterValues;
    OpenCLArray<cl_float>* sumBuffer;
1108
    OpenCLArray<cl_float>* energy;
1109
1110
    OpenCLArray<mm_float4>* uniformRandoms;
    OpenCLArray<mm_int4>* randomSeed;
1111
1112
1113
    OpenCLParameterSet* perDofValues;
    mutable std::vector<std::vector<cl_float> > localPerDofValues;
    std::vector<std::vector<cl::Kernel> > kernels;
1114
    cl::Kernel sumEnergyKernel, randomKernel;
1115
1116
1117
1118
    std::vector<CustomIntegrator::ComputationType> stepType;
    std::vector<bool> needsForces;
    std::vector<bool> needsEnergy;
    std::vector<bool> invalidatesForces;
1119
    std::vector<bool> merged;
1120
    std::vector<int> forceGroup;
1121
1122
    std::vector<int> requiredGaussian;
    std::vector<int> requiredUniform;
1123
    std::vector<std::string> parameterNames;
1124
1125
};

1126
1127
1128
1129
1130
1131
/**
 * 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),
1132
            hasInitializedKernels(false), atomGroups(NULL) {
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
    }
    ~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;
1152
    OpenCLArray<cl_int>* atomGroups;
1153
    cl::Kernel kernel;
1154
1155
1156
1157
1158
1159
1160
1161
};

/**
 * 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),
1162
            hasInitializedKernels(false), savedPositions(NULL), moleculeAtoms(NULL), moleculeStartIndex(NULL) {
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
    }
    ~OpenCLApplyMonteCarloBarostatKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param barostat   the MonteCarloBarostat this kernel will be used for
     */
    void initialize(const System& system, const MonteCarloBarostat& barostat);
    /**
     * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
     * 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 scale      the scale factor by which to multiply particle positions
     */
    void scaleCoordinates(ContextImpl& context, double scale);
    /**
     * 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;
    OpenCLArray<mm_float4>* savedPositions;
    OpenCLArray<cl_int>* moleculeAtoms;
    OpenCLArray<cl_int>* moleculeStartIndex;
    cl::Kernel kernel;
1197
};
1198
1199
1200
1201
1202
1203

/**
 * This kernel is invoked to calculate the kinetic energy of the system.
 */
class OpenCLCalcKineticEnergyKernel : public CalcKineticEnergyKernel {
public:
1204
    OpenCLCalcKineticEnergyKernel(std::string name, const Platform& platform, OpenCLContext& cl) : CalcKineticEnergyKernel(name, platform), cl(cl) {
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
    }
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     */
    void initialize(const System& system);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
     */
    double execute(ContextImpl& context);
private:
1219
    OpenCLContext& cl;
1220
1221
1222
    std::vector<double> masses;
};

1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
/**
 * 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;
    OpenCLArray<mm_float4>* cmMomentum;
    cl::Kernel kernel1, kernel2;
};
1250
1251
1252
1253

} // namespace OpenMM

#endif /*OPENMM_OPENCLKERNELS_H_*/