ReferenceKernels.h 45.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#ifndef OPENMM_REFERENCEKERNELS_H_
#define OPENMM_REFERENCEKERNELS_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-2009 Stanford University and the Authors.      *
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
 * Authors: Peter Eastman                                                     *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

35
#include "ReferencePlatform.h"
36
#include "openmm/kernels.h"
37
#include "SimTKUtilities/SimTKOpenMMRealType.h"
38
#include "SimTKReference/ReferenceNeighborList.h"
39
#include "lepton/ExpressionProgram.h"
40

41
class CpuObc;
Mark Friedrichs's avatar
Mark Friedrichs committed
42
class CpuGBVI;
43
class ReferenceAndersenThermostat;
44
class ReferenceCustomCompoundBondIxn;
45
class ReferenceCustomHbondIxn;
46
class ReferenceBrownianDynamics;
47
class ReferenceStochasticDynamics;
48
class ReferenceConstraintAlgorithm;
49
class ReferenceMonteCarloBarostat;
50
class ReferenceVariableStochasticDynamics;
51
class ReferenceVariableVerletDynamics;
52
class ReferenceVerletDynamics;
53
class ReferenceCustomDynamics;
54

55
56
namespace OpenMM {

57
/**
58
59
60
 * 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.
61
 */
62
class ReferenceCalcForcesAndEnergyKernel : public CalcForcesAndEnergyKernel {
63
public:
64
    ReferenceCalcForcesAndEnergyKernel(std::string name, const Platform& platform) : CalcForcesAndEnergyKernel(name, platform) {
65
66
67
68
69
70
71
72
    }
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     */
    void initialize(const System& system);
    /**
73
     * This is called at the beginning of each force/energy computation, before calcForcesAndEnergy() has been called on
74
75
     * any ForceImpl.
     *
76
77
78
     * @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
79
     * @param groups        a set of bit flags for which force groups to include
80
     */
81
    void beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
82
    /**
83
     * This is called at the end of each force/energy computation, after calcForcesAndEnergy() has been called on
84
85
     * every ForceImpl.
     *
86
87
88
     * @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
89
     * @param groups        a set of bit flags for which force groups to include
90
     * @return the potential energy of the system.  This value is added to all values returned by ForceImpls'
91
     * calcForcesAndEnergy() methods.  That is, each force kernel may <i>either</i> return its contribution to the
92
     * energy directly, <i>or</i> add it to an internal buffer so that it will be included here.
93
     */
94
    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
95
96
private:
    std::vector<RealVec> savedForces;
97
98
};

99
/**
100
101
 * This kernel provides methods for setting and retrieving various state data: time, positions,
 * velocities, and forces.
102
 */
103
class ReferenceUpdateStateDataKernel : public UpdateStateDataKernel {
104
public:
105
    ReferenceUpdateStateDataKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : UpdateStateDataKernel(name, platform), data(data) {
106
107
108
109
110
111
112
113
114
115
116
117
    }
    /**
     * 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
     */
118
    double getTime(const ContextImpl& context) const;
119
120
121
122
123
    /**
     * Set the current time (in picoseconds).
     *
     * @param context    the context in which to execute this kernel
     */
124
    void setTime(ContextImpl& context, double time);
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
    /**
     * 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);
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
    /**
     * 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
171
172
173
174
175
176
177
178
179
180
181
182
    /**
     * 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);
183
184
185
186
private:
    ReferencePlatform::PlatformData& data;
};

187
188
189
190
191
192
/**
 * This kernel modifies the positions of particles to enforce distance constraints.
 */
class ReferenceApplyConstraintsKernel : public ApplyConstraintsKernel {
public:
    ReferenceApplyConstraintsKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) :
193
            ApplyConstraintsKernel(name, platform), data(data), constraints(0), constraintDistances(0), constraintIndices(0) {
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
    }
    ~ReferenceApplyConstraintsKernel();
    /**
     * 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:
    ReferencePlatform::PlatformData& data;
    ReferenceConstraintAlgorithm* constraints;
212
213
    std::vector<RealOpenMM> masses;
    std::vector<RealOpenMM> inverseMasses;
214
215
216
217
218
    RealOpenMM* constraintDistances;
    int** constraintIndices;
    int numConstraints;
};

219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
/**
 * This kernel recomputes the positions of virtual sites.
 */
class ReferenceVirtualSitesKernel : public VirtualSitesKernel {
public:
    ReferenceVirtualSitesKernel(std::string name, const Platform& platform) : VirtualSitesKernel(name, platform) {
    }
    /**
     * 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);
};

240
/**
241
 * This kernel is invoked by HarmonicBondForce to calculate the forces acting on the system and the energy of the system.
242
 */
243
class ReferenceCalcHarmonicBondForceKernel : public CalcHarmonicBondForceKernel {
244
public:
245
    ReferenceCalcHarmonicBondForceKernel(std::string name, const Platform& platform) : CalcHarmonicBondForceKernel(name, platform) {
246
    }
247
    ~ReferenceCalcHarmonicBondForceKernel();
248
    /**
249
     * Initialize the kernel.
250
     * 
251
     * @param system     the System this kernel will be applied to
252
253
254
255
     * @param force      the HarmonicBondForce this kernel will be used for
     */
    void initialize(const System& system, const HarmonicBondForce& force);
    /**
256
257
258
259
260
261
     * 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
262
     */
263
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
264
265
266
267
268
269
private:
    int numBonds;
    int **bondIndexArray;
    RealOpenMM **bondParamArray;
};

270
271
272
273
274
275
276
277
278
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 ReferenceCalcCustomBondForceKernel : public CalcCustomBondForceKernel {
public:
    ReferenceCalcCustomBondForceKernel(std::string name, const Platform& platform) : CalcCustomBondForceKernel(name, platform) {
    }
    ~ReferenceCalcCustomBondForceKernel();
    /**
     * 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);
    /**
286
     * Execute the kernel to calculate the forces and/or energy.
287
     *
288
289
290
291
     * @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
292
     */
293
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
294
295
296
297
298
299
300
301
private:
    int numBonds;
    int **bondIndexArray;
    RealOpenMM **bondParamArray;
    Lepton::ExpressionProgram energyExpression, forceExpression;
    std::vector<std::string> parameterNames, globalParameterNames;
};

302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
/**
 * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
public:
    ReferenceCalcHarmonicAngleForceKernel(std::string name, const Platform& platform) : CalcHarmonicAngleForceKernel(name, platform) {
    }
    ~ReferenceCalcHarmonicAngleForceKernel();
    /**
     * 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);
    /**
318
319
320
321
322
323
     * 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
324
     */
325
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
326
327
328
329
330
331
private:
    int numAngles;
    int **angleIndexArray;
    RealOpenMM **angleParamArray;
};

332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
/**
 * This kernel is invoked by CustomAngleForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcCustomAngleForceKernel : public CalcCustomAngleForceKernel {
public:
    ReferenceCalcCustomAngleForceKernel(std::string name, const Platform& platform) : CalcCustomAngleForceKernel(name, platform) {
    }
    ~ReferenceCalcCustomAngleForceKernel();
    /**
     * 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);
    /**
348
     * Execute the kernel to calculate the forces and/or energy.
349
     *
350
351
352
353
     * @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
354
     */
355
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
356
357
358
359
360
361
362
363
private:
    int numAngles;
    int **angleIndexArray;
    RealOpenMM **angleParamArray;
    Lepton::ExpressionProgram energyExpression, forceExpression;
    std::vector<std::string> parameterNames, globalParameterNames;
};

364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
/**
 * This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
public:
    ReferenceCalcPeriodicTorsionForceKernel(std::string name, const Platform& platform) : CalcPeriodicTorsionForceKernel(name, platform) {
    }
    ~ReferenceCalcPeriodicTorsionForceKernel();
    /**
     * 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);
    /**
380
381
382
383
384
385
     * 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
386
     */
387
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
private:
    int numTorsions;
    int **torsionIndexArray;
    RealOpenMM **torsionParamArray;
};

/**
 * This kernel is invoked by RBTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
public:
    ReferenceCalcRBTorsionForceKernel(std::string name, const Platform& platform) : CalcRBTorsionForceKernel(name, platform) {
    }
    ~ReferenceCalcRBTorsionForceKernel();
    /**
     * 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);
    /**
410
411
412
413
414
415
     * 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
416
     */
417
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
418
419
420
421
422
423
private:
    int numTorsions;
    int **torsionIndexArray;
    RealOpenMM **torsionParamArray;
};

424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
/**
 * This kernel is invoked by CMAPTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcCMAPTorsionForceKernel : public CalcCMAPTorsionForceKernel {
public:
    ReferenceCalcCMAPTorsionForceKernel(std::string name, const Platform& platform) : CalcCMAPTorsionForceKernel(name, platform) {
    }
    /**
     * 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);
    /**
439
     * Execute the kernel to calculate the forces and/or energy.
440
     *
441
442
443
444
     * @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
445
     */
446
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
447
448
449
450
451
452
private:
    std::vector<std::vector<std::vector<RealOpenMM> > > coeff;
    std::vector<int> torsionMaps;
    std::vector<std::vector<int> > torsionIndices;
};

453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
/**
 * This kernel is invoked by CustomTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcCustomTorsionForceKernel : public CalcCustomTorsionForceKernel {
public:
    ReferenceCalcCustomTorsionForceKernel(std::string name, const Platform& platform) : CalcCustomTorsionForceKernel(name, platform) {
    }
    ~ReferenceCalcCustomTorsionForceKernel();
    /**
     * 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);
    /**
469
     * Execute the kernel to calculate the forces and/or energy.
470
     *
471
472
473
474
     * @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
475
     */
476
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
477
478
479
480
481
482
483
484
private:
    int numTorsions;
    int **torsionIndexArray;
    RealOpenMM **torsionParamArray;
    Lepton::ExpressionProgram energyExpression, forceExpression;
    std::vector<std::string> parameterNames, globalParameterNames;
};

485
486
487
488
489
490
491
492
493
494
495
496
497
/**
 * This kernel is invoked by NonbondedForce to calculate the forces acting on the system.
 */
class ReferenceCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
public:
    ReferenceCalcNonbondedForceKernel(std::string name, const Platform& platform) : CalcNonbondedForceKernel(name, platform) {
    }
    ~ReferenceCalcNonbondedForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param force      the NonbondedForce this kernel will be used for
498
     */
499
    void initialize(const System& system, const NonbondedForce& force);
500
    /**
501
502
503
504
505
     * 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
506
     * @param includeReciprocal  true if reciprocal space interactions should be included
507
     * @return the potential energy due to the force
508
     */
509
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal);
510
511
512
513
514
515
516
    /**
     * 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);
517
private:
Peter Eastman's avatar
Peter Eastman committed
518
    int numParticles, num14;
519
    int **exclusionArray, **bonded14IndexArray;
Peter Eastman's avatar
Peter Eastman committed
520
    RealOpenMM **particleParamArray, **bonded14ParamArray;
521
    RealOpenMM nonbondedCutoff, rfDielectric, ewaldAlpha, dispersionCoefficient;
522
    int kmax[3], gridSize[3];
523
524
525
    std::vector<std::set<int> > exclusions;
    NonbondedMethod nonbondedMethod;
    NeighborList* neighborList;
526
527
};

528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
/**
 * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system.
 */
class ReferenceCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
public:
    ReferenceCalcCustomNonbondedForceKernel(std::string name, const Platform& platform) : CalcCustomNonbondedForceKernel(name, platform) {
    }
    ~ReferenceCalcCustomNonbondedForceKernel();
    /**
     * 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);
    /**
544
     * Execute the kernel to calculate the forces and/or energy.
545
     *
546
547
548
549
     * @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
550
     */
551
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
552
private:
553
554
555
    int numParticles;
    int **exclusionArray;
    RealOpenMM **particleParamArray;
556
557
558
    RealOpenMM nonbondedCutoff, periodicBoxSize[3];
    std::vector<std::set<int> > exclusions;
    Lepton::ExpressionProgram energyExpression, forceExpression;
559
    std::vector<std::string> parameterNames, globalParameterNames;
560
561
562
563
    NonbondedMethod nonbondedMethod;
    NeighborList* neighborList;
};

564
/**
565
 * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
566
 */
567
class ReferenceCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
568
public:
569
    ReferenceCalcGBSAOBCForceKernel(std::string name, const Platform& platform) : CalcGBSAOBCForceKernel(name, platform) {
570
    }
571
    ~ReferenceCalcGBSAOBCForceKernel();
572
    /**
573
     * Initialize the kernel.
574
     * 
575
     * @param system     the System this kernel will be applied to
576
     * @param force      the GBSAOBCForce this kernel will be used for
577
     */
578
    void initialize(const System& system, const GBSAOBCForce& force);
579
    /**
580
581
582
583
584
585
     * 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
586
     */
587
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
588
589
590
private:
    CpuObc* obc;
    std::vector<RealOpenMM> charges;
591
    bool isPeriodic;
592
593
};

Mark Friedrichs's avatar
Mark Friedrichs committed
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
/**
 * This kernel is invoked by GBVIForce to calculate the forces acting on the system.
 */
class ReferenceCalcGBVIForceKernel : public CalcGBVIForceKernel {
public:
    ReferenceCalcGBVIForceKernel(std::string name, const Platform& platform) : CalcGBVIForceKernel(name, platform) {
    }
    ~ReferenceCalcGBVIForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system       the System this kernel will be applied to
     * @param force        the GBVIForce this kernel will be used for
     * @param scaled radii the scaled radii (Eq. 5 of Labute paper)
     */
    void initialize(const System& system, const GBVIForce& force, const std::vector<double> & scaledRadii);
    /**
611
612
613
614
615
616
     * 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
Mark Friedrichs's avatar
Mark Friedrichs committed
617
     */
618
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
Mark Friedrichs's avatar
Mark Friedrichs committed
619
620
621
private:
    CpuGBVI * gbvi;
    std::vector<RealOpenMM> charges;
622
    bool isPeriodic;
Mark Friedrichs's avatar
Mark Friedrichs committed
623
624
};

625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
/**
 * This kernel is invoked by CustomGBForce to calculate the forces acting on the system.
 */
class ReferenceCalcCustomGBForceKernel : public CalcCustomGBForceKernel {
public:
    ReferenceCalcCustomGBForceKernel(std::string name, const Platform& platform) : CalcCustomGBForceKernel(name, platform) {
    }
    ~ReferenceCalcCustomGBForceKernel();
    /**
     * 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);
    /**
641
     * Execute the kernel to calculate the forces and/or energy.
642
     *
643
644
645
646
     * @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
647
     */
648
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
649
650
private:
    int numParticles;
651
    bool isPeriodic;
652
    RealOpenMM **particleParamArray;
653
    RealOpenMM nonbondedCutoff;
654
655
656
    std::vector<std::set<int> > exclusions;
    std::vector<std::string> particleParameterNames, globalParameterNames, valueNames;
    std::vector<Lepton::ExpressionProgram> valueExpressions;
657
    std::vector<std::vector<Lepton::ExpressionProgram> > valueDerivExpressions;
658
    std::vector<std::vector<Lepton::ExpressionProgram> > valueGradientExpressions;
659
660
661
    std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
    std::vector<Lepton::ExpressionProgram> energyExpressions;
    std::vector<std::vector<Lepton::ExpressionProgram> > energyDerivExpressions;
662
    std::vector<std::vector<Lepton::ExpressionProgram> > energyGradientExpressions;
663
664
665
666
667
    std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
    NonbondedMethod nonbondedMethod;
    NeighborList* neighborList;
};

668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
/**
 * This kernel is invoked by CustomExternalForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcCustomExternalForceKernel : public CalcCustomExternalForceKernel {
public:
    ReferenceCalcCustomExternalForceKernel(std::string name, const Platform& platform) : CalcCustomExternalForceKernel(name, platform) {
    }
    ~ReferenceCalcCustomExternalForceKernel();
    /**
     * 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);
    /**
684
     * Execute the kernel to calculate the forces and/or energy.
685
     *
686
687
688
689
     * @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
690
     */
691
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
692
693
694
695
696
697
698
699
private:
    int numParticles;
    std::vector<int> particles;
    RealOpenMM **particleParamArray;
    Lepton::ExpressionProgram energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ;
    std::vector<std::string> parameterNames, globalParameterNames;
};

700
701
702
703
704
/**
 * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system.
 */
class ReferenceCalcCustomHbondForceKernel : public CalcCustomHbondForceKernel {
public:
705
    ReferenceCalcCustomHbondForceKernel(std::string name, const Platform& platform) : CalcCustomHbondForceKernel(name, platform), ixn(NULL) {
706
707
708
709
710
711
712
713
714
715
    }
    ~ReferenceCalcCustomHbondForceKernel();
    /**
     * 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);
    /**
716
     * Execute the kernel to calculate the forces and/or energy.
717
     *
718
719
720
721
     * @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
722
     */
723
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
724
725
private:
    int numDonors, numAcceptors, numParticles;
726
    bool isPeriodic;
727
728
    int **exclusionArray;
    RealOpenMM **donorParamArray, **acceptorParamArray;
729
    RealOpenMM nonbondedCutoff;
730
    ReferenceCustomHbondIxn* ixn;
731
    std::vector<std::set<int> > exclusions;
732
    std::vector<std::string> globalParameterNames;
733
734
};

735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
/**
 * This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system.
 */
class ReferenceCalcCustomCompoundBondForceKernel : public CalcCustomCompoundBondForceKernel {
public:
    ReferenceCalcCustomCompoundBondForceKernel(std::string name, const Platform& platform) : CalcCustomCompoundBondForceKernel(name, platform), ixn(NULL) {
    }
    ~ReferenceCalcCustomCompoundBondForceKernel();
    /**
     * 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, numParticles;
    RealOpenMM **bondParamArray;
    ReferenceCustomCompoundBondIxn* ixn;
    std::vector<std::string> globalParameterNames;
};

766
767
768
769
770
/**
 * This kernel is invoked by VerletIntegrator to take one time step.
 */
class ReferenceIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
public:
771
    ReferenceIntegrateVerletStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVerletStepKernel(name, platform),
772
        data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
773
    }
774
    ~ReferenceIntegrateVerletStepKernel();
775
    /**
776
     * Initialize the kernel.
777
     * 
778
779
     * @param system     the System this kernel will be applied to
     * @param integrator the VerletIntegrator this kernel will be used for
780
     */
781
    void initialize(const System& system, const VerletIntegrator& integrator);
782
783
784
    /**
     * Execute the kernel.
     * 
785
786
     * @param context    the context in which to execute this kernel
     * @param integrator the VerletIntegrator this kernel is being used for
787
     */
788
    void execute(ContextImpl& context, const VerletIntegrator& integrator);
789
private:
790
    ReferencePlatform::PlatformData& data;
791
    ReferenceVerletDynamics* dynamics;
792
    ReferenceConstraintAlgorithm* constraints;
793
    std::vector<RealOpenMM> masses;
794
    RealOpenMM* constraintDistances;
795
796
797
    int** constraintIndices;
    int numConstraints;
    double prevStepSize;
798
799
800
801
802
803
804
};

/**
 * This kernel is invoked by LangevinIntegrator to take one time step.
 */
class ReferenceIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
public:
805
    ReferenceIntegrateLangevinStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateLangevinStepKernel(name, platform),
806
        data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
807
    }
808
    ~ReferenceIntegrateLangevinStepKernel();
809
    /**
Peter Eastman's avatar
Peter Eastman committed
810
     * Initialize the kernel, setting up the particle masses.
811
     * 
812
813
     * @param system     the System this kernel will be applied to
     * @param integrator the LangevinIntegrator this kernel will be used for
814
     */
815
    void initialize(const System& system, const LangevinIntegrator& integrator);
816
817
818
    /**
     * Execute the kernel.
     * 
819
820
     * @param context    the context in which to execute this kernel
     * @param integrator the LangevinIntegrator this kernel is being used for
821
     */
822
    void execute(ContextImpl& context, const LangevinIntegrator& integrator);
823
private:
824
    ReferencePlatform::PlatformData& data;
825
    ReferenceStochasticDynamics* dynamics;
826
    ReferenceConstraintAlgorithm* constraints;
827
    std::vector<RealOpenMM> masses;
828
    RealOpenMM* constraintDistances;
829
830
831
    int** constraintIndices;
    int numConstraints;
    double prevTemp, prevFriction, prevStepSize;
832
833
834
835
836
837
838
};

/**
 * This kernel is invoked by BrownianIntegrator to take one time step.
 */
class ReferenceIntegrateBrownianStepKernel : public IntegrateBrownianStepKernel {
public:
839
    ReferenceIntegrateBrownianStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateBrownianStepKernel(name, platform),
840
        data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
841
    }
842
    ~ReferenceIntegrateBrownianStepKernel();
843
    /**
844
     * Initialize the kernel.
845
     * 
846
847
     * @param system     the System this kernel will be applied to
     * @param integrator the BrownianIntegrator this kernel will be used for
848
     */
849
    void initialize(const System& system, const BrownianIntegrator& integrator);
850
851
852
    /**
     * Execute the kernel.
     * 
853
854
     * @param context    the context in which to execute this kernel
     * @param integrator the BrownianIntegrator this kernel is being used for
855
     */
856
    void execute(ContextImpl& context, const BrownianIntegrator& integrator);
857
private:
858
    ReferencePlatform::PlatformData& data;
859
    ReferenceBrownianDynamics* dynamics;
860
    ReferenceConstraintAlgorithm* constraints;
861
    std::vector<RealOpenMM> masses;
862
    RealOpenMM* constraintDistances;
863
864
865
    int** constraintIndices;
    int numConstraints;
    double prevTemp, prevFriction, prevStepSize;
866
867
};

868
869
870
871
872
873
/**
 * This kernel is invoked by VariableLangevinIntegrator to take one time step.
 */
class ReferenceIntegrateVariableLangevinStepKernel : public IntegrateVariableLangevinStepKernel {
public:
    ReferenceIntegrateVariableLangevinStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVariableLangevinStepKernel(name, platform),
874
        data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
    }
    ~ReferenceIntegrateVariableLangevinStepKernel();
    /**
     * Initialize the kernel.
     *
     * @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 VariableLangevinIntegrator& integrator);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
     * @param integrator the LangevinIntegrator this kernel is being used for
     * @param maxTime    the maximum time beyond which the simulation should not be advanced
890
     * @return the size of the step that was taken
891
     */
892
    double execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
893
894
895
896
private:
    ReferencePlatform::PlatformData& data;
    ReferenceVariableStochasticDynamics* dynamics;
    ReferenceConstraintAlgorithm* constraints;
897
    std::vector<RealOpenMM> masses;
898
899
900
901
902
903
    RealOpenMM* constraintDistances;
    int** constraintIndices;
    int numConstraints;
    double prevTemp, prevFriction, prevErrorTol;
};

904
905
906
907
908
909
/**
 * This kernel is invoked by VariableVerletIntegrator to take one time step.
 */
class ReferenceIntegrateVariableVerletStepKernel : public IntegrateVariableVerletStepKernel {
public:
    ReferenceIntegrateVariableVerletStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVariableVerletStepKernel(name, platform),
910
        data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
911
912
913
914
915
916
917
918
919
920
921
922
923
924
    }
    ~ReferenceIntegrateVariableVerletStepKernel();
    /**
     * 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
925
     * @param maxTime    the maximum time beyond which the simulation should not be advanced
926
     * @return the size of the step that was taken
927
     */
928
    double execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
929
930
931
932
private:
    ReferencePlatform::PlatformData& data;
    ReferenceVariableVerletDynamics* dynamics;
    ReferenceConstraintAlgorithm* constraints;
933
    std::vector<RealOpenMM> masses;
934
935
936
    RealOpenMM* constraintDistances;
    int** constraintIndices;
    int numConstraints;
937
    double prevErrorTol;
938
939
};

940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
/**
 * This kernel is invoked by CustomIntegrator to take one time step.
 */
class ReferenceIntegrateCustomStepKernel : public IntegrateCustomStepKernel {
public:
    ReferenceIntegrateCustomStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateCustomStepKernel(name, platform),
        data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
    }
    ~ReferenceIntegrateCustomStepKernel();
    /**
     * 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:
    ReferencePlatform::PlatformData& data;
    ReferenceCustomDynamics* dynamics;
    ReferenceConstraintAlgorithm* constraints;
    std::vector<RealOpenMM> masses, globalValues;
    std::vector<std::vector<OpenMM::RealVec> > perDofValues; 
    RealOpenMM* constraintDistances;
    int** constraintIndices;
    int numConstraints;
};

1008
/**
Peter Eastman's avatar
Peter Eastman committed
1009
 * This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
1010
1011
1012
 */
class ReferenceApplyAndersenThermostatKernel : public ApplyAndersenThermostatKernel {
public:
1013
    ReferenceApplyAndersenThermostatKernel(std::string name, const Platform& platform) : ApplyAndersenThermostatKernel(name, platform), thermostat(0) {
1014
    }
1015
    ~ReferenceApplyAndersenThermostatKernel();
1016
    /**
1017
     * Initialize the kernel.
1018
     * 
1019
1020
     * @param system     the System this kernel will be applied to
     * @param thermostat the AndersenThermostat this kernel will be used for
1021
     */
1022
    void initialize(const System& system, const AndersenThermostat& thermostat);
1023
1024
1025
    /**
     * Execute the kernel.
     * 
1026
     * @param context    the context in which to execute this kernel
1027
     */
1028
    void execute(ContextImpl& context);
1029
1030
private:
    ReferenceAndersenThermostat* thermostat;
1031
    std::vector<std::vector<int> > particleGroups;
1032
    std::vector<RealOpenMM> masses;
1033
1034
};

1035
1036
1037
1038
1039
/**
 * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
 */
class ReferenceApplyMonteCarloBarostatKernel : public ApplyMonteCarloBarostatKernel {
public:
1040
    ReferenceApplyMonteCarloBarostatKernel(std::string name, const Platform& platform) : ApplyMonteCarloBarostatKernel(name, platform), barostat(NULL) {
1041
1042
1043
1044
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
    }
    ~ReferenceApplyMonteCarloBarostatKernel();
    /**
     * 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:
    ReferenceMonteCarloBarostat* barostat;
};

1071
1072
1073
1074
1075
/**
 * This kernel is invoked to calculate the kinetic energy of the system.
 */
class ReferenceCalcKineticEnergyKernel : public CalcKineticEnergyKernel {
public:
1076
    ReferenceCalcKineticEnergyKernel(std::string name, const Platform& platform) : CalcKineticEnergyKernel(name, platform) {
1077
1078
    }
    /**
1079
     * Initialize the kernel.
1080
     * 
1081
     * @param system     the System this kernel will be applied to
1082
     */
1083
    void initialize(const System& system);
1084
1085
1086
    /**
     * Execute the kernel.
     * 
1087
     * @param context    the context in which to execute this kernel
1088
     */
1089
    double execute(ContextImpl& context);
1090
1091
private:
    std::vector<double> masses;
1092
1093
};

1094
1095
1096
1097
1098
/**
 * This kernel is invoked to remove center of mass motion from the system.
 */
class ReferenceRemoveCMMotionKernel : public RemoveCMMotionKernel {
public:
1099
    ReferenceRemoveCMMotionKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : RemoveCMMotionKernel(name, platform), data(data) {
1100
1101
    }
    /**
Peter Eastman's avatar
Peter Eastman committed
1102
     * Initialize the kernel, setting up the particle masses.
1103
     * 
1104
1105
     * @param system     the System this kernel will be applied to
     * @param force      the CMMotionRemover this kernel will be used for
1106
     */
1107
    void initialize(const System& system, const CMMotionRemover& force);
1108
1109
1110
    /**
     * Execute the kernel.
     * 
1111
     * @param context    the context in which to execute this kernel
1112
     */
1113
    void execute(ContextImpl& context);
1114
private:
1115
    ReferencePlatform::PlatformData& data;
1116
    std::vector<double> masses;
1117
    int frequency;
1118
1119
};

1120
1121
1122
} // namespace OpenMM

#endif /*OPENMM_REFERENCEKERNELS_H_*/