OpenCLKernels.h 66.9 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-2017 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
#include "openmm/kernels.h"
37
38
39
#include "openmm/internal/CompiledExpressionSet.h"
#include "openmm/internal/CustomIntegratorUtilities.h"
#include "lepton/CompiledExpression.h"
40
41
42
43
#include "openmm/System.h"

namespace OpenMM {

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

/**
89
90
 * This kernel provides methods for setting and retrieving various state data: time, positions,
 * velocities, and forces.
91
 */
92
class OpenCLUpdateStateDataKernel : public UpdateStateDataKernel {
93
public:
94
    OpenCLUpdateStateDataKernel(std::string name, const Platform& platform, OpenCLContext& cl) : UpdateStateDataKernel(name, platform), cl(cl) {
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
    }
    /**
     * 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);
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
141
142
143
    /**
     * 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);
144
145
146
147
148
149
    /**
     * Get the current derivatives of the energy with respect to context parameters.
     *
     * @param derivs  on exit, this contains the derivatives
     */
    void getEnergyParameterDerivatives(ContextImpl& context, std::map<std::string, double>& derivs);
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
    /**
     * 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
     */
165
    void setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c);
Peter Eastman's avatar
Peter Eastman committed
166
167
168
169
170
171
172
173
174
175
176
177
    /**
     * 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);
178
private:
179
    OpenCLContext& cl;
180
};
181

182
183
184
185
186
/**
 * This kernel modifies the positions of particles to enforce distance constraints.
 */
class OpenCLApplyConstraintsKernel : public ApplyConstraintsKernel {
public:
187
188
    OpenCLApplyConstraintsKernel(std::string name, const Platform& platform, OpenCLContext& cl) : ApplyConstraintsKernel(name, platform),
            cl(cl), hasInitializedKernel(false) {
189
190
191
192
193
194
195
196
197
198
199
200
201
202
    }
    /**
     * 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);
203
204
205
206
207
208
209
    /**
     * 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);
210
211
private:
    OpenCLContext& cl;
212
213
    bool hasInitializedKernel;
    cl::Kernel applyDeltasKernel;
214
215
};

216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
/**
 * 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;
};

239
240
241
242
243
/**
 * 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:
244
    OpenCLCalcHarmonicBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcHarmonicBondForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
245
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
246
247
248
249
250
251
252
253
254
255
    }
    ~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);
    /**
256
     * Execute the kernel to calculate the forces and/or energy.
257
     *
258
259
260
261
     * @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
262
     */
263
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
264
265
266
267
268
269
270
    /**
     * 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);
271
private:
272
    class ForceInfo;
273
    int numBonds;
274
    bool hasInitializedKernel;
275
    OpenCLContext& cl;
276
    ForceInfo* info;
277
    const System& system;
278
    OpenCLArray* params;
279
280
};

281
282
283
284
285
/**
 * 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:
286
    OpenCLCalcCustomBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomBondForceKernel(name, platform),
287
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
288
289
290
291
292
293
294
295
296
297
    }
    ~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);
    /**
298
     * Execute the kernel to calculate the forces and/or energy.
299
     *
300
301
302
303
     * @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
304
     */
305
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
306
307
308
309
310
311
312
    /**
     * 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);
313
private:
314
    class ForceInfo;
315
316
317
    int numBonds;
    bool hasInitializedKernel;
    OpenCLContext& cl;
318
    ForceInfo* info;
319
    const System& system;
320
    OpenCLParameterSet* params;
321
    OpenCLArray* globals;
322
323
324
325
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

326
327
328
329
330
/**
 * 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:
331
    OpenCLCalcHarmonicAngleForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcHarmonicAngleForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
332
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
333
334
335
336
337
338
339
340
341
342
    }
    ~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);
    /**
343
     * Execute the kernel to calculate the forces and/or energy.
344
     *
345
346
347
348
     * @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
349
     */
350
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
351
352
353
354
355
356
357
    /**
     * 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);
358
private:
359
    class ForceInfo;
360
    int numAngles;
361
    bool hasInitializedKernel;
362
    OpenCLContext& cl;
363
    ForceInfo* info;
364
    const System& system;
365
    OpenCLArray* params;
366
367
};

368
369
370
371
372
/**
 * 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:
373
    OpenCLCalcCustomAngleForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomAngleForceKernel(name, platform),
374
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
375
376
377
378
379
380
381
382
383
384
    }
    ~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);
    /**
385
     * Execute the kernel to calculate the forces and/or energy.
386
     *
387
388
389
390
     * @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
391
     */
392
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
393
394
395
396
397
398
399
    /**
     * 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);
400
private:
401
    class ForceInfo;
402
403
404
    int numAngles;
    bool hasInitializedKernel;
    OpenCLContext& cl;
405
    ForceInfo* info;
406
    const System& system;
407
    OpenCLParameterSet* params;
408
    OpenCLArray* globals;
409
410
411
412
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

413
414
415
416
417
/**
 * 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:
418
    OpenCLCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcPeriodicTorsionForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
419
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
420
421
422
423
424
425
426
427
428
429
    }
    ~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);
    /**
430
     * Execute the kernel to calculate the forces and/or energy.
431
     *
432
433
434
435
     * @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
436
     */
437
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
438
439
440
441
442
443
444
    /**
     * 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);
445
private:
446
    class ForceInfo;
447
    int numTorsions;
448
    bool hasInitializedKernel;
449
    OpenCLContext& cl;
450
    ForceInfo* info;
451
    const System& system;
452
    OpenCLArray* params;
453
454
};

455
456
457
458
459
/**
 * 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:
460
    OpenCLCalcRBTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcRBTorsionForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
461
            hasInitializedKernel(false), cl(cl), system(system), params(NULL) {
462
463
464
465
466
467
468
469
470
471
    }
    ~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);
    /**
472
     * Execute the kernel to calculate the forces and/or energy.
473
     *
474
475
476
477
     * @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
478
     */
479
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
480
481
482
483
484
485
486
    /**
     * 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);
487
private:
488
    class ForceInfo;
489
    int numTorsions;
490
    bool hasInitializedKernel;
491
    OpenCLContext& cl;
492
    ForceInfo* info;
493
    const System& system;
494
    OpenCLArray* params;
495
496
};

497
498
499
500
501
/**
 * 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:
502
    OpenCLCalcCMAPTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCMAPTorsionForceKernel(name, platform),
503
            hasInitializedKernel(false), cl(cl), system(system), coefficients(NULL), mapPositions(NULL), torsionMaps(NULL) {
504
505
506
507
508
509
510
511
512
513
    }
    ~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);
    /**
514
     * Execute the kernel to calculate the forces and/or energy.
515
     *
516
517
518
519
     * @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
520
     */
521
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
522
523
524
525
526
527
528
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CMAPTorsionForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force);
529
private:
530
    class ForceInfo;
531
532
533
    int numTorsions;
    bool hasInitializedKernel;
    OpenCLContext& cl;
534
    ForceInfo* info;
535
    const System& system;
536
    std::vector<mm_int2> mapPositionsVec;
537
538
539
    OpenCLArray* coefficients;
    OpenCLArray* mapPositions;
    OpenCLArray* torsionMaps;
540
541
};

542
543
544
545
546
/**
 * 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:
547
    OpenCLCalcCustomTorsionForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomTorsionForceKernel(name, platform),
548
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
549
550
551
552
553
554
555
556
557
558
    }
    ~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);
    /**
559
     * Execute the kernel to calculate the forces and/or energy.
560
     *
561
562
563
564
     * @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
565
     */
566
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
567
568
569
570
571
572
573
    /**
     * 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);
574
private:
575
    class ForceInfo;
576
577
578
    int numTorsions;
    bool hasInitializedKernel;
    OpenCLContext& cl;
579
    ForceInfo* info;
580
    const System& system;
581
    OpenCLParameterSet* params;
582
    OpenCLArray* globals;
583
584
585
586
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

587
588
589
590
591
/**
 * This kernel is invoked by NonbondedForce to calculate the forces acting on the system.
 */
class OpenCLCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
public:
592
    OpenCLCalcNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcNonbondedForceKernel(name, platform),
Peter Eastman's avatar
Peter Eastman committed
593
            hasInitializedKernel(false), cl(cl), sigmaEpsilon(NULL), exceptionParams(NULL), cosSinSums(NULL), pmeGrid(NULL),
594
            pmeGrid2(NULL), pmeBsplineModuliX(NULL), pmeBsplineModuliY(NULL), pmeBsplineModuliZ(NULL), pmeBsplineTheta(NULL),
595
            pmeAtomRange(NULL), pmeAtomGridIndex(NULL), pmeEnergyBuffer(NULL), sort(NULL), fft(NULL), pmeio(NULL) {
596
597
598
599
600
601
602
603
604
605
    }
    ~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);
    /**
606
     * Execute the kernel to calculate the forces and/or energy.
607
     *
608
609
610
     * @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
611
612
     * @param includeDirect  true if direct space interactions should be included
     * @param includeReciprocal  true if reciprocal space interactions should be included
613
     * @return the potential energy due to the force
614
     */
615
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal);
616
617
618
619
620
621
622
    /**
     * 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);
623
624
625
626
627
628
629
630
631
    /**
     * Get the parameters being used for PME.
     * 
     * @param alpha   the separation parameter
     * @param nx      the number of grid points along the X axis
     * @param ny      the number of grid points along the Y axis
     * @param nz      the number of grid points along the Z axis
     */
    void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
632
private:
633
634
635
636
637
638
639
640
641
    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";}
642
    };
643
    class ForceInfo;
644
645
646
    class PmeIO;
    class PmePreComputation;
    class PmePostComputation;
647
648
    class SyncQueuePreComputation;
    class SyncQueuePostComputation;
649
    OpenCLContext& cl;
650
    ForceInfo* info;
651
    bool hasInitializedKernel;
652
653
654
655
656
657
658
659
660
661
662
    OpenCLArray* sigmaEpsilon;
    OpenCLArray* exceptionParams;
    OpenCLArray* cosSinSums;
    OpenCLArray* pmeGrid;
    OpenCLArray* pmeGrid2;
    OpenCLArray* pmeBsplineModuliX;
    OpenCLArray* pmeBsplineModuliY;
    OpenCLArray* pmeBsplineModuliZ;
    OpenCLArray* pmeBsplineTheta;
    OpenCLArray* pmeAtomRange;
    OpenCLArray* pmeAtomGridIndex;
663
    OpenCLArray* pmeEnergyBuffer;
664
    OpenCLSort* sort;
665
666
    cl::CommandQueue pmeQueue;
    cl::Event pmeSyncEvent;
667
    OpenCLFFT3D* fft;
668
669
    Kernel cpuPme;
    PmeIO* pmeio;
670
    SyncQueuePostComputation* syncQueue;
671
672
    cl::Kernel ewaldSumsKernel;
    cl::Kernel ewaldForcesKernel;
673
674
    cl::Kernel pmeGridIndexKernel;
    cl::Kernel pmeAtomRangeKernel;
675
    cl::Kernel pmeZIndexKernel;
676
677
    cl::Kernel pmeUpdateBsplinesKernel;
    cl::Kernel pmeSpreadChargeKernel;
678
    cl::Kernel pmeFinishSpreadChargeKernel;
679
    cl::Kernel pmeConvolutionKernel;
680
    cl::Kernel pmeEvalEnergyKernel;
681
682
    cl::Kernel pmeInterpolateForceKernel;
    std::map<std::string, std::string> pmeDefines;
683
    std::vector<std::pair<int, int> > exceptionAtoms;
684
    double ewaldSelfEnergy, dispersionCoefficient, alpha;
685
    int gridSizeX, gridSizeY, gridSizeZ;
686
    bool hasCoulomb, hasLJ, usePmeQueue;
687
    NonbondedMethod nonbondedMethod;
688
    static const int PmeOrder = 5;
689
690
};

691
692
693
694
695
/**
 * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
public:
696
    OpenCLCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomNonbondedForceKernel(name, platform),
697
            cl(cl), params(NULL), globals(NULL), interactionGroupData(NULL), forceCopy(NULL), system(system), hasInitializedKernel(false) {
698
699
700
701
702
703
704
705
706
707
    }
    ~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);
    /**
708
     * Execute the kernel to calculate the forces and/or energy.
709
     *
710
711
712
713
     * @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
714
     */
715
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
716
717
718
719
720
721
722
    /**
     * 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);
723
private:
724
    class ForceInfo;
725
    void initInteractionGroups(const CustomNonbondedForce& force, const std::string& interactionSource, const std::vector<std::string>& tableTypes);
726
    OpenCLContext& cl;
727
    ForceInfo* info;
728
    OpenCLParameterSet* params;
729
    OpenCLArray* globals;
730
731
732
    OpenCLArray* interactionGroupData;
    cl::Kernel interactionGroupKernel;
    std::vector<void*> interactionGroupArgs;
733
734
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
735
    std::vector<OpenCLArray*> tabulatedFunctions;
736
    double longRangeCoefficient;
737
    std::vector<double> longRangeCoefficientDerivs;
738
739
    bool hasInitializedLongRangeCorrection, hasInitializedKernel;
    int numGroupThreadBlocks;
740
    CustomNonbondedForce* forceCopy;
741
    const System& system;
742
};
743
744
745
746
747
748

/**
 * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
 */
class OpenCLCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
public:
749
    OpenCLCalcGBSAOBCForceKernel(std::string name, const Platform& platform, OpenCLContext& cl) : CalcGBSAOBCForceKernel(name, platform), cl(cl),
750
751
            hasCreatedKernels(false), params(NULL), bornSum(NULL), longBornSum(NULL), bornRadii(NULL), bornForce(NULL),
            longBornForce(NULL), obcChain(NULL) {
752
753
754
755
756
757
758
759
760
761
    }
    ~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);
    /**
762
     * Execute the kernel to calculate the forces and/or energy.
763
     *
764
765
766
767
     * @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
768
     */
769
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
770
771
772
773
774
775
776
    /**
     * 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);
777
private:
778
    class ForceInfo;
779
    double prefactor, surfaceAreaFactor, cutoff;
780
    bool hasCreatedKernels;
781
    int maxTiles;
782
    OpenCLContext& cl;
783
    ForceInfo* info;
784
785
786
787
788
789
790
    OpenCLArray* params;
    OpenCLArray* bornSum;
    OpenCLArray* longBornSum;
    OpenCLArray* bornRadii;
    OpenCLArray* bornForce;
    OpenCLArray* longBornForce;
    OpenCLArray* obcChain;
791
792
    cl::Kernel computeBornSumKernel;
    cl::Kernel reduceBornSumKernel;
793
794
    cl::Kernel force1Kernel;
    cl::Kernel reduceBornForceKernel;
795
};
796

797
798
799
800
801
/**
 * This kernel is invoked by CustomGBForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomGBForceKernel : public CalcCustomGBForceKernel {
public:
802
    OpenCLCalcCustomGBForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomGBForceKernel(name, platform),
803
            hasInitializedKernels(false), cl(cl), params(NULL), computedValues(NULL), energyDerivs(NULL), energyDerivChain(NULL), longEnergyDerivs(NULL), globals(NULL),
804
            valueBuffers(NULL), longValueBuffers(NULL), system(system) {
805
806
807
808
809
810
811
812
813
814
    }
    ~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);
    /**
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 CustomGBForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomGBForce& force);
830
private:
831
    class ForceInfo;
832
    double cutoff;
833
    bool hasInitializedKernels, needParameterGradient, needEnergyParamDerivs;
834
    int maxTiles, numComputedValues;
835
    OpenCLContext& cl;
836
    ForceInfo* info;
837
838
    OpenCLParameterSet* params;
    OpenCLParameterSet* computedValues;
839
    OpenCLParameterSet* energyDerivs;
840
    OpenCLParameterSet* energyDerivChain;
841
842
    std::vector<OpenCLParameterSet*> dValuedParam;
    std::vector<OpenCLArray*> dValue0dParam;
843
844
845
846
    OpenCLArray* longEnergyDerivs;
    OpenCLArray* globals;
    OpenCLArray* valueBuffers;
    OpenCLArray* longValueBuffers;
847
848
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
849
    std::vector<OpenCLArray*> tabulatedFunctions;
Peter Eastman's avatar
Peter Eastman committed
850
    std::vector<bool> pairValueUsesParam, pairEnergyUsesParam, pairEnergyUsesValue;
851
    const System& system;
Peter Eastman's avatar
Peter Eastman committed
852
    cl::Kernel pairValueKernel, perParticleValueKernel, pairEnergyKernel, perParticleEnergyKernel, gradientChainRuleKernel;
853
854
    std::string pairValueSrc, pairEnergySrc;
    std::map<std::string, std::string> pairValueDefines, pairEnergyDefines;
855
856
};

857
858
859
860
861
/**
 * 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:
862
    OpenCLCalcCustomExternalForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomExternalForceKernel(name, platform),
863
            hasInitializedKernel(false), cl(cl), system(system), params(NULL), globals(NULL) {
864
865
866
867
868
869
870
871
872
873
    }
    ~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);
    /**
874
     * Execute the kernel to calculate the forces and/or energy.
875
     *
876
877
878
879
     * @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
880
     */
881
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
882
883
884
885
886
887
888
    /**
     * 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);
889
private:
890
    class ForceInfo;
891
892
893
    int numParticles;
    bool hasInitializedKernel;
    OpenCLContext& cl;
894
    ForceInfo* info;
895
    const System& system;
896
    OpenCLParameterSet* params;
897
    OpenCLArray* globals;
898
899
900
901
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
};

902
903
904
905
906
/**
 * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomHbondForceKernel : public CalcCustomHbondForceKernel {
public:
907
    OpenCLCalcCustomHbondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomHbondForceKernel(name, platform),
908
            hasInitializedKernel(false), cl(cl), donorParams(NULL), acceptorParams(NULL), donors(NULL), acceptors(NULL),
909
            donorBufferIndices(NULL), acceptorBufferIndices(NULL), globals(NULL), donorExclusions(NULL), acceptorExclusions(NULL), system(system) {
910
911
912
913
914
915
916
917
918
919
    }
    ~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);
    /**
920
     * Execute the kernel to calculate the forces and/or energy.
921
     *
922
923
924
925
     * @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
926
     */
927
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
928
929
930
931
932
933
934
    /**
     * 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);
935
private:
936
    class ForceInfo;
937
938
939
    int numDonors, numAcceptors;
    bool hasInitializedKernel;
    OpenCLContext& cl;
940
    ForceInfo* info;
941
942
    OpenCLParameterSet* donorParams;
    OpenCLParameterSet* acceptorParams;
943
944
945
946
947
948
949
    OpenCLArray* globals;
    OpenCLArray* donors;
    OpenCLArray* acceptors;
    OpenCLArray* donorBufferIndices;
    OpenCLArray* acceptorBufferIndices;
    OpenCLArray* donorExclusions;
    OpenCLArray* acceptorExclusions;
950
951
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
952
    std::vector<OpenCLArray*> tabulatedFunctions;
953
    const System& system;
954
    cl::Kernel donorKernel, acceptorKernel;
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
/**
 * This kernel is invoked by CustomCentroidBondForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomCentroidBondForceKernel : public CalcCustomCentroidBondForceKernel {
public:
    OpenCLCalcCustomCentroidBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomCentroidBondForceKernel(name, platform),
            cl(cl), params(NULL), globals(NULL), groupParticles(NULL), groupWeights(NULL), groupOffsets(NULL), groupForces(NULL), bondGroups(NULL), centerPositions(NULL), system(system) {
    }
    ~OpenCLCalcCustomCentroidBondForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomCentroidBondForce this kernel will be used for
     */
    void initialize(const System& system, const CustomCentroidBondForce& 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 CustomCentroidBondForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force);

private:
991
    class ForceInfo;
992
    int numGroups, numBonds;
993
    bool needEnergyParamDerivs;
994
    OpenCLContext& cl;
995
    ForceInfo* info;
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
    OpenCLParameterSet* params;
    OpenCLArray* globals;
    OpenCLArray* groupParticles;
    OpenCLArray* groupWeights;
    OpenCLArray* groupOffsets;
    OpenCLArray* groupForces;
    OpenCLArray* bondGroups;
    OpenCLArray* centerPositions;
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
    std::vector<OpenCLArray*> tabulatedFunctions;
    cl::Kernel computeCentersKernel, groupForcesKernel, applyForcesKernel;
    const System& system;
};

1011
1012
1013
1014
1015
/**
 * This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system.
 */
class OpenCLCalcCustomCompoundBondForceKernel : public CalcCustomCompoundBondForceKernel {
public:
1016
    OpenCLCalcCustomCompoundBondForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomCompoundBondForceKernel(name, platform),
1017
            cl(cl), params(NULL), globals(NULL), system(system) {
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
    }
    ~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);
1036
1037
1038
1039
1040
1041
1042
1043
    /**
     * 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);

1044
private:
1045
    class ForceInfo;
1046
1047
    int numBonds;
    OpenCLContext& cl;
1048
    ForceInfo* info;
1049
    OpenCLParameterSet* params;
1050
    OpenCLArray* globals;
1051
1052
    std::vector<std::string> globalParamNames;
    std::vector<cl_float> globalParamValues;
1053
    std::vector<OpenCLArray*> tabulatedFunctions;
1054
    const System& system;
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
/**
 * 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:
1093
    class ForceInfo;
1094
    OpenCLContext& cl;
1095
    ForceInfo* info;
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
    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;
};

1120
1121
1122
1123
1124
1125
/**
 * This kernel is invoked by GayBerneForce to calculate the forces acting on the system.
 */
class OpenCLCalcGayBerneForceKernel : public CalcGayBerneForceKernel {
public:
    OpenCLCalcGayBerneForceKernel(std::string name, const Platform& platform, OpenCLContext& cl) : CalcGayBerneForceKernel(name, platform), cl(cl),
1126
1127
1128
1129
            hasInitializedKernels(false), sortedParticles(NULL), axisParticleIndices(NULL), sigParams(NULL), epsParams(NULL), scale(NULL), exceptionParticles(NULL),
            exceptionParams(NULL), aMatrix(NULL),
            bMatrix(NULL), gMatrix(NULL), exclusions(NULL), exclusionStartIndex(NULL), blockCenter(NULL), blockBoundingBox(NULL), neighbors(NULL),
            neighborIndex(NULL), neighborBlockCount(NULL), sortedPos(NULL), torque(NULL) {
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
    }
    ~OpenCLCalcGayBerneForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the GayBerneForce this kernel will be used for
     */
    void initialize(const System& system, const GayBerneForce& 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
     * @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 GayBerneForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const GayBerneForce& force);
private:
1155
    class ForceInfo;
1156
1157
1158
    class ReorderListener;
    void sortAtoms();
    OpenCLContext& cl;
1159
    ForceInfo* info;
1160
    bool hasInitializedKernels;
1161
    int numRealParticles, maxNeighborBlocks;
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
    GayBerneForce::NonbondedMethod nonbondedMethod;
    OpenCLArray* sortedParticles;
    OpenCLArray* axisParticleIndices;
    OpenCLArray* sigParams;
    OpenCLArray* epsParams;
    OpenCLArray* scale;
    OpenCLArray* exceptionParticles;
    OpenCLArray* exceptionParams;
    OpenCLArray* aMatrix;
    OpenCLArray* bMatrix;
    OpenCLArray* gMatrix;
    OpenCLArray* exclusions;
    OpenCLArray* exclusionStartIndex;
    OpenCLArray* blockCenter;
    OpenCLArray* blockBoundingBox;
1177
1178
1179
    OpenCLArray* neighbors;
    OpenCLArray* neighborIndex;
    OpenCLArray* neighborBlockCount;
1180
1181
1182
1183
1184
    OpenCLArray* sortedPos;
    OpenCLArray* torque;
    std::vector<bool> isRealParticle;
    std::vector<std::pair<int, int> > exceptionAtoms;
    std::vector<std::pair<int, int> > excludedPairs;
1185
    cl::Kernel framesKernel, blockBoundsKernel, neighborsKernel, forceKernel, torqueKernel;
1186
1187
};

1188
1189
1190
1191
1192
/**
 * This kernel is invoked by VerletIntegrator to take one time step.
 */
class OpenCLIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
public:
1193
    OpenCLIntegrateVerletStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateVerletStepKernel(name, platform), cl(cl),
1194
            hasInitializedKernels(false) {
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
    }
    ~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);
1211
1212
1213
1214
1215
1216
1217
    /**
     * 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);
1218
private:
1219
    OpenCLContext& cl;
1220
    bool hasInitializedKernels;
1221
1222
1223
1224
1225
1226
1227
1228
1229
    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),
1230
            hasInitializedKernels(false), params(NULL) {
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
    }
    ~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);
1247
1248
1249
1250
1251
1252
1253
    /**
     * 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);
1254
1255
1256
private:
    OpenCLContext& cl;
    double prevTemp, prevFriction, prevStepSize;
1257
    bool hasInitializedKernels;
1258
    OpenCLArray* params;
1259
    cl::Kernel kernel1, kernel2;
1260
1261
};

1262
1263
1264
1265
1266
/**
 * This kernel is invoked by BrownianIntegrator to take one time step.
 */
class OpenCLIntegrateBrownianStepKernel : public IntegrateBrownianStepKernel {
public:
1267
1268
    OpenCLIntegrateBrownianStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateBrownianStepKernel(name, platform), cl(cl),
            hasInitializedKernels(false), prevTemp(-1), prevFriction(-1), prevStepSize(-1) {
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
    }
    ~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);
1285
1286
1287
1288
1289
1290
1291
    /**
     * 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);
1292
1293
1294
1295
1296
1297
private:
    OpenCLContext& cl;
    double prevTemp, prevFriction, prevStepSize;
    bool hasInitializedKernels;
    cl::Kernel kernel1, kernel2;
};
1298
1299
1300
1301
1302
1303
1304

/**
 * 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),
1305
            hasInitializedKernels(false) {
1306
1307
1308
1309
1310
1311
    }
    ~OpenCLIntegrateVariableVerletStepKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
1312
     * @param integrator the VariableVerletIntegrator this kernel will be used for
1313
1314
1315
1316
1317
1318
     */
    void initialize(const System& system, const VariableVerletIntegrator& integrator);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
1319
     * @param integrator the VariableVerletIntegrator this kernel is being used for
1320
     * @param maxTime    the maximum time beyond which the simulation should not be advanced
1321
     * @return the size of the step that was taken
1322
     */
1323
    double execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
1324
1325
1326
1327
1328
1329
1330
    /**
     * 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);
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
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),
1344
            hasInitializedKernels(false), params(NULL) {
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
    }
    ~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
1360
     * @return the size of the step that was taken
1361
     */
1362
    double execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
1363
1364
1365
1366
1367
1368
1369
    /**
     * 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);
1370
1371
1372
1373
private:
    OpenCLContext& cl;
    bool hasInitializedKernels;
    int blockSize;
1374
    OpenCLArray* params;
1375
    cl::Kernel kernel1, kernel2, selectSizeKernel;
1376
1377
    double prevTemp, prevFriction, prevErrorTol;
};
1378

1379
1380
1381
1382
1383
/**
 * This kernel is invoked by CustomIntegrator to take one time step.
 */
class OpenCLIntegrateCustomStepKernel : public IntegrateCustomStepKernel {
public:
1384
    enum GlobalTargetType {DT, VARIABLE, PARAMETER};
1385
    OpenCLIntegrateCustomStepKernel(std::string name, const Platform& platform, OpenCLContext& cl) : IntegrateCustomStepKernel(name, platform), cl(cl),
1386
            hasInitializedKernels(false), localValuesAreCurrent(false), globalValues(NULL), sumBuffer(NULL), summedValue(NULL), uniformRandoms(NULL),
1387
            randomSeed(NULL), perDofEnergyParamDerivs(NULL), perDofValues(NULL), needsEnergyParamDerivs(false) {
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
    }
    ~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);
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
    /**
     * 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);
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
    /**
     * 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:
1450
    class ReorderListener;
1451
    class GlobalTarget;
1452
    class DerivFunction;
1453
    std::string createPerDofComputation(const std::string& variable, const Lepton::ParsedExpression& expr, int component, CustomIntegrator& integrator, const std::string& forceName, const std::string& energyName);
1454
    void prepareForComputation(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid);
1455
1456
    Lepton::ExpressionTreeNode replaceDerivFunctions(const Lepton::ExpressionTreeNode& node, OpenMM::ContextImpl& context);
    void findExpressionsForDerivs(const Lepton::ExpressionTreeNode& node, std::vector<std::pair<Lepton::ExpressionTreeNode, std::string> >& variableNodes);
1457
    void recordGlobalValue(double value, GlobalTarget target, CustomIntegrator& integrator);
1458
    void recordChangedParameters(ContextImpl& context);
1459
    bool evaluateCondition(int step);
1460
    OpenCLContext& cl;
1461
    double energy;
1462
    float energyFloat;
1463
    int numGlobalVariables;
1464
    bool hasInitializedKernels, deviceValuesAreCurrent, deviceGlobalsAreCurrent, modifiesParameters, keNeedsForce, hasAnyConstraints, needsEnergyParamDerivs;
1465
    mutable bool localValuesAreCurrent;
1466
1467
    OpenCLArray* globalValues;
    OpenCLArray* sumBuffer;
1468
    OpenCLArray* summedValue;
1469
1470
    OpenCLArray* uniformRandoms;
    OpenCLArray* randomSeed;
1471
    OpenCLArray* perDofEnergyParamDerivs;
1472
1473
    std::map<int, OpenCLArray*> savedForces;
    std::set<int> validSavedForces;
1474
    OpenCLParameterSet* perDofValues;
1475
1476
    mutable std::vector<std::vector<cl_float> > localPerDofValuesFloat;
    mutable std::vector<std::vector<cl_double> > localPerDofValuesDouble;
1477
1478
1479
1480
    std::map<std::string, double> energyParamDerivs;
    std::vector<std::string> perDofEnergyParamDerivNames;
    std::vector<cl_float> localPerDofEnergyParamDerivsFloat;
    std::vector<cl_double> localPerDofEnergyParamDerivsDouble;
1481
1482
1483
    std::vector<float> globalValuesFloat;
    std::vector<double> globalValuesDouble;
    std::vector<double> initialGlobalVariables;
1484
    std::vector<std::vector<cl::Kernel> > kernels;
Peter Eastman's avatar
Peter Eastman committed
1485
    cl::Kernel randomKernel, kineticEnergyKernel, sumKineticEnergyKernel;
1486
    std::vector<CustomIntegrator::ComputationType> stepType;
1487
1488
1489
1490
    std::vector<CustomIntegratorUtilities::Comparison> comparisons;
    std::vector<std::vector<Lepton::CompiledExpression> > globalExpressions;
    CompiledExpressionSet expressionSet;
    std::vector<bool> needsGlobals;
1491
1492
    std::vector<bool> needsForces;
    std::vector<bool> needsEnergy;
1493
    std::vector<bool> computeBothForceAndEnergy;
1494
    std::vector<bool> invalidatesForces;
1495
    std::vector<bool> merged;
1496
1497
    std::vector<int> forceGroupFlags;
    std::vector<int> blockEnd;
1498
1499
    std::vector<int> requiredGaussian;
    std::vector<int> requiredUniform;
1500
1501
1502
1503
    std::vector<int> stepEnergyVariableIndex;
    std::vector<int> globalVariableIndex;
    std::vector<int> parameterVariableIndex;
    int gaussianVariableIndex, uniformVariableIndex, dtVariableIndex;
1504
    std::vector<std::string> parameterNames;
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
    std::vector<GlobalTarget> stepTarget;
};

class OpenCLIntegrateCustomStepKernel::GlobalTarget {
public:
    OpenCLIntegrateCustomStepKernel::GlobalTargetType type;
    int variableIndex;
    GlobalTarget() {
    }
    GlobalTarget(OpenCLIntegrateCustomStepKernel::GlobalTargetType type, int variableIndex) : type(type), variableIndex(variableIndex) {
    }
1516
1517
};

1518
1519
1520
1521
1522
1523
/**
 * 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),
1524
            hasInitializedKernels(false), atomGroups(NULL) {
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
    }
    ~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;
1544
    OpenCLArray* atomGroups;
1545
    cl::Kernel kernel;
1546
1547
1548
1549
1550
1551
1552
1553
};

/**
 * 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),
1554
            hasInitializedKernels(false), savedPositions(NULL), moleculeAtoms(NULL), moleculeStartIndex(NULL) {
1555
1556
1557
1558
1559
1560
1561
1562
    }
    ~OpenCLApplyMonteCarloBarostatKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param barostat   the MonteCarloBarostat this kernel will be used for
     */
1563
    void initialize(const System& system, const Force& barostat);
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
    /**
     * 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
     */
1576
    void scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ);
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
    /**
     * 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;
1588
1589
1590
    OpenCLArray* savedPositions;
    OpenCLArray* moleculeAtoms;
    OpenCLArray* moleculeStartIndex;
1591
    cl::Kernel kernel;
1592
    std::vector<int> lastAtomOrder;
1593
};
1594

1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
/**
 * 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;
1619
    OpenCLArray* cmMomentum;
1620
1621
    cl::Kernel kernel1, kernel2;
};
1622
1623
1624
1625

} // namespace OpenMM

#endif /*OPENMM_OPENCLKERNELS_H_*/