ReferenceKernels.h 44.7 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;
171
172
173
174
private:
    ReferencePlatform::PlatformData& data;
};

175
176
177
178
179
180
/**
 * 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) :
181
            ApplyConstraintsKernel(name, platform), data(data), constraints(0), constraintDistances(0), constraintIndices(0) {
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
    }
    ~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;
200
201
    std::vector<RealOpenMM> masses;
    std::vector<RealOpenMM> inverseMasses;
202
203
204
205
206
    RealOpenMM* constraintDistances;
    int** constraintIndices;
    int numConstraints;
};

207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
/**
 * 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);
};

228
/**
229
 * This kernel is invoked by HarmonicBondForce to calculate the forces acting on the system and the energy of the system.
230
 */
231
class ReferenceCalcHarmonicBondForceKernel : public CalcHarmonicBondForceKernel {
232
public:
233
    ReferenceCalcHarmonicBondForceKernel(std::string name, const Platform& platform) : CalcHarmonicBondForceKernel(name, platform) {
234
    }
235
    ~ReferenceCalcHarmonicBondForceKernel();
236
    /**
237
     * Initialize the kernel.
238
     * 
239
     * @param system     the System this kernel will be applied to
240
241
242
243
     * @param force      the HarmonicBondForce this kernel will be used for
     */
    void initialize(const System& system, const HarmonicBondForce& force);
    /**
244
245
246
247
248
249
     * 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
250
     */
251
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
252
253
254
255
256
257
private:
    int numBonds;
    int **bondIndexArray;
    RealOpenMM **bondParamArray;
};

258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
/**
 * 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);
    /**
274
     * Execute the kernel to calculate the forces and/or energy.
275
     *
276
277
278
279
     * @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
280
     */
281
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
282
283
284
285
286
287
288
289
private:
    int numBonds;
    int **bondIndexArray;
    RealOpenMM **bondParamArray;
    Lepton::ExpressionProgram energyExpression, forceExpression;
    std::vector<std::string> parameterNames, globalParameterNames;
};

290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
/**
 * 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);
    /**
306
307
308
309
310
311
     * 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
312
     */
313
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
314
315
316
317
318
319
private:
    int numAngles;
    int **angleIndexArray;
    RealOpenMM **angleParamArray;
};

320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
/**
 * 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);
    /**
336
     * Execute the kernel to calculate the forces and/or energy.
337
     *
338
339
340
341
     * @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
342
     */
343
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
344
345
346
347
348
349
350
351
private:
    int numAngles;
    int **angleIndexArray;
    RealOpenMM **angleParamArray;
    Lepton::ExpressionProgram energyExpression, forceExpression;
    std::vector<std::string> parameterNames, globalParameterNames;
};

352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
/**
 * 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);
    /**
368
369
370
371
372
373
     * 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
374
     */
375
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
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);
    /**
398
399
400
401
402
403
     * 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
404
     */
405
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
406
407
408
409
410
411
private:
    int numTorsions;
    int **torsionIndexArray;
    RealOpenMM **torsionParamArray;
};

412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
/**
 * 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);
    /**
427
     * Execute the kernel to calculate the forces and/or energy.
428
     *
429
430
431
432
     * @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
433
     */
434
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
435
436
437
438
439
440
private:
    std::vector<std::vector<std::vector<RealOpenMM> > > coeff;
    std::vector<int> torsionMaps;
    std::vector<std::vector<int> > torsionIndices;
};

441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
/**
 * 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);
    /**
457
     * Execute the kernel to calculate the forces and/or energy.
458
     *
459
460
461
462
     * @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
463
     */
464
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
465
466
467
468
469
470
471
472
private:
    int numTorsions;
    int **torsionIndexArray;
    RealOpenMM **torsionParamArray;
    Lepton::ExpressionProgram energyExpression, forceExpression;
    std::vector<std::string> parameterNames, globalParameterNames;
};

473
474
475
476
477
478
479
480
481
482
483
484
485
/**
 * 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
486
     */
487
    void initialize(const System& system, const NonbondedForce& force);
488
    /**
489
490
491
492
493
     * 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
494
     * @param includeReciprocal  true if reciprocal space interactions should be included
495
     * @return the potential energy due to the force
496
     */
497
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal);
498
private:
Peter Eastman's avatar
Peter Eastman committed
499
    int numParticles, num14;
500
    int **exclusionArray, **bonded14IndexArray;
Peter Eastman's avatar
Peter Eastman committed
501
    RealOpenMM **particleParamArray, **bonded14ParamArray;
502
    RealOpenMM nonbondedCutoff, rfDielectric, ewaldAlpha, dispersionCoefficient;
503
    int kmax[3], gridSize[3];
504
505
506
    std::vector<std::set<int> > exclusions;
    NonbondedMethod nonbondedMethod;
    NeighborList* neighborList;
507
508
};

509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
/**
 * 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);
    /**
525
     * Execute the kernel to calculate the forces and/or energy.
526
     *
527
528
529
530
     * @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
531
     */
532
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
533
private:
534
535
536
    int numParticles;
    int **exclusionArray;
    RealOpenMM **particleParamArray;
537
538
539
    RealOpenMM nonbondedCutoff, periodicBoxSize[3];
    std::vector<std::set<int> > exclusions;
    Lepton::ExpressionProgram energyExpression, forceExpression;
540
    std::vector<std::string> parameterNames, globalParameterNames;
541
542
543
544
    NonbondedMethod nonbondedMethod;
    NeighborList* neighborList;
};

545
/**
546
 * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
547
 */
548
class ReferenceCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
549
public:
550
    ReferenceCalcGBSAOBCForceKernel(std::string name, const Platform& platform) : CalcGBSAOBCForceKernel(name, platform) {
551
    }
552
    ~ReferenceCalcGBSAOBCForceKernel();
553
    /**
554
     * Initialize the kernel.
555
     * 
556
     * @param system     the System this kernel will be applied to
557
     * @param force      the GBSAOBCForce this kernel will be used for
558
     */
559
    void initialize(const System& system, const GBSAOBCForce& force);
560
    /**
561
562
563
564
565
566
     * 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
567
     */
568
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
569
570
571
private:
    CpuObc* obc;
    std::vector<RealOpenMM> charges;
572
    bool isPeriodic;
573
574
};

Mark Friedrichs's avatar
Mark Friedrichs committed
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
/**
 * 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);
    /**
592
593
594
595
596
597
     * 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
598
     */
599
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
Mark Friedrichs's avatar
Mark Friedrichs committed
600
601
602
private:
    CpuGBVI * gbvi;
    std::vector<RealOpenMM> charges;
603
    bool isPeriodic;
Mark Friedrichs's avatar
Mark Friedrichs committed
604
605
};

606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
/**
 * 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);
    /**
622
     * Execute the kernel to calculate the forces and/or energy.
623
     *
624
625
626
627
     * @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
628
     */
629
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
630
631
private:
    int numParticles;
632
    bool isPeriodic;
633
    RealOpenMM **particleParamArray;
634
    RealOpenMM nonbondedCutoff;
635
636
637
    std::vector<std::set<int> > exclusions;
    std::vector<std::string> particleParameterNames, globalParameterNames, valueNames;
    std::vector<Lepton::ExpressionProgram> valueExpressions;
638
    std::vector<std::vector<Lepton::ExpressionProgram> > valueDerivExpressions;
639
    std::vector<std::vector<Lepton::ExpressionProgram> > valueGradientExpressions;
640
641
642
    std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
    std::vector<Lepton::ExpressionProgram> energyExpressions;
    std::vector<std::vector<Lepton::ExpressionProgram> > energyDerivExpressions;
643
    std::vector<std::vector<Lepton::ExpressionProgram> > energyGradientExpressions;
644
645
646
647
648
    std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
    NonbondedMethod nonbondedMethod;
    NeighborList* neighborList;
};

649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
/**
 * 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);
    /**
665
     * Execute the kernel to calculate the forces and/or energy.
666
     *
667
668
669
670
     * @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
671
     */
672
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
673
674
675
676
677
678
679
680
private:
    int numParticles;
    std::vector<int> particles;
    RealOpenMM **particleParamArray;
    Lepton::ExpressionProgram energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ;
    std::vector<std::string> parameterNames, globalParameterNames;
};

681
682
683
684
685
/**
 * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system.
 */
class ReferenceCalcCustomHbondForceKernel : public CalcCustomHbondForceKernel {
public:
686
    ReferenceCalcCustomHbondForceKernel(std::string name, const Platform& platform) : CalcCustomHbondForceKernel(name, platform), ixn(NULL) {
687
688
689
690
691
692
693
694
695
696
    }
    ~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);
    /**
697
     * Execute the kernel to calculate the forces and/or energy.
698
     *
699
700
701
702
     * @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
703
     */
704
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
705
706
private:
    int numDonors, numAcceptors, numParticles;
707
    bool isPeriodic;
708
709
    int **exclusionArray;
    RealOpenMM **donorParamArray, **acceptorParamArray;
710
    RealOpenMM nonbondedCutoff;
711
    ReferenceCustomHbondIxn* ixn;
712
    std::vector<std::set<int> > exclusions;
713
    std::vector<std::string> globalParameterNames;
714
715
};

716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
/**
 * 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;
};

747
748
749
750
751
/**
 * This kernel is invoked by VerletIntegrator to take one time step.
 */
class ReferenceIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
public:
752
    ReferenceIntegrateVerletStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVerletStepKernel(name, platform),
753
        data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
754
    }
755
    ~ReferenceIntegrateVerletStepKernel();
756
    /**
757
     * Initialize the kernel.
758
     * 
759
760
     * @param system     the System this kernel will be applied to
     * @param integrator the VerletIntegrator this kernel will be used for
761
     */
762
    void initialize(const System& system, const VerletIntegrator& integrator);
763
764
765
    /**
     * Execute the kernel.
     * 
766
767
     * @param context    the context in which to execute this kernel
     * @param integrator the VerletIntegrator this kernel is being used for
768
     */
769
    void execute(ContextImpl& context, const VerletIntegrator& integrator);
770
private:
771
    ReferencePlatform::PlatformData& data;
772
    ReferenceVerletDynamics* dynamics;
773
    ReferenceConstraintAlgorithm* constraints;
774
    std::vector<RealOpenMM> masses;
775
    RealOpenMM* constraintDistances;
776
777
778
    int** constraintIndices;
    int numConstraints;
    double prevStepSize;
779
780
781
782
783
784
785
};

/**
 * This kernel is invoked by LangevinIntegrator to take one time step.
 */
class ReferenceIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
public:
786
    ReferenceIntegrateLangevinStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateLangevinStepKernel(name, platform),
787
        data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
788
    }
789
    ~ReferenceIntegrateLangevinStepKernel();
790
    /**
Peter Eastman's avatar
Peter Eastman committed
791
     * Initialize the kernel, setting up the particle masses.
792
     * 
793
794
     * @param system     the System this kernel will be applied to
     * @param integrator the LangevinIntegrator this kernel will be used for
795
     */
796
    void initialize(const System& system, const LangevinIntegrator& integrator);
797
798
799
    /**
     * Execute the kernel.
     * 
800
801
     * @param context    the context in which to execute this kernel
     * @param integrator the LangevinIntegrator this kernel is being used for
802
     */
803
    void execute(ContextImpl& context, const LangevinIntegrator& integrator);
804
private:
805
    ReferencePlatform::PlatformData& data;
806
    ReferenceStochasticDynamics* dynamics;
807
    ReferenceConstraintAlgorithm* constraints;
808
    std::vector<RealOpenMM> masses;
809
    RealOpenMM* constraintDistances;
810
811
812
    int** constraintIndices;
    int numConstraints;
    double prevTemp, prevFriction, prevStepSize;
813
814
815
816
817
818
819
};

/**
 * This kernel is invoked by BrownianIntegrator to take one time step.
 */
class ReferenceIntegrateBrownianStepKernel : public IntegrateBrownianStepKernel {
public:
820
    ReferenceIntegrateBrownianStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateBrownianStepKernel(name, platform),
821
        data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
822
    }
823
    ~ReferenceIntegrateBrownianStepKernel();
824
    /**
825
     * Initialize the kernel.
826
     * 
827
828
     * @param system     the System this kernel will be applied to
     * @param integrator the BrownianIntegrator this kernel will be used for
829
     */
830
    void initialize(const System& system, const BrownianIntegrator& integrator);
831
832
833
    /**
     * Execute the kernel.
     * 
834
835
     * @param context    the context in which to execute this kernel
     * @param integrator the BrownianIntegrator this kernel is being used for
836
     */
837
    void execute(ContextImpl& context, const BrownianIntegrator& integrator);
838
private:
839
    ReferencePlatform::PlatformData& data;
840
    ReferenceBrownianDynamics* dynamics;
841
    ReferenceConstraintAlgorithm* constraints;
842
    std::vector<RealOpenMM> masses;
843
    RealOpenMM* constraintDistances;
844
845
846
    int** constraintIndices;
    int numConstraints;
    double prevTemp, prevFriction, prevStepSize;
847
848
};

849
850
851
852
853
854
/**
 * 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),
855
        data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
    }
    ~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
871
     * @return the size of the step that was taken
872
     */
873
    double execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
874
875
876
877
private:
    ReferencePlatform::PlatformData& data;
    ReferenceVariableStochasticDynamics* dynamics;
    ReferenceConstraintAlgorithm* constraints;
878
    std::vector<RealOpenMM> masses;
879
880
881
882
883
884
    RealOpenMM* constraintDistances;
    int** constraintIndices;
    int numConstraints;
    double prevTemp, prevFriction, prevErrorTol;
};

885
886
887
888
889
890
/**
 * 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),
891
        data(data), dynamics(0), constraints(0), constraintDistances(0), constraintIndices(0) {
892
893
894
895
896
897
898
899
900
901
902
903
904
905
    }
    ~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
906
     * @param maxTime    the maximum time beyond which the simulation should not be advanced
907
     * @return the size of the step that was taken
908
     */
909
    double execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
910
911
912
913
private:
    ReferencePlatform::PlatformData& data;
    ReferenceVariableVerletDynamics* dynamics;
    ReferenceConstraintAlgorithm* constraints;
914
    std::vector<RealOpenMM> masses;
915
916
917
    RealOpenMM* constraintDistances;
    int** constraintIndices;
    int numConstraints;
918
    double prevErrorTol;
919
920
};

921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
/**
 * 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;
};

989
/**
Peter Eastman's avatar
Peter Eastman committed
990
 * This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
991
992
993
 */
class ReferenceApplyAndersenThermostatKernel : public ApplyAndersenThermostatKernel {
public:
994
    ReferenceApplyAndersenThermostatKernel(std::string name, const Platform& platform) : ApplyAndersenThermostatKernel(name, platform), thermostat(0) {
995
    }
996
    ~ReferenceApplyAndersenThermostatKernel();
997
    /**
998
     * Initialize the kernel.
999
     * 
1000
1001
     * @param system     the System this kernel will be applied to
     * @param thermostat the AndersenThermostat this kernel will be used for
1002
     */
1003
    void initialize(const System& system, const AndersenThermostat& thermostat);
1004
1005
1006
    /**
     * Execute the kernel.
     * 
1007
     * @param context    the context in which to execute this kernel
1008
     */
1009
    void execute(ContextImpl& context);
1010
1011
private:
    ReferenceAndersenThermostat* thermostat;
1012
    std::vector<std::vector<int> > particleGroups;
1013
    std::vector<RealOpenMM> masses;
1014
1015
};

1016
1017
1018
1019
1020
/**
 * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
 */
class ReferenceApplyMonteCarloBarostatKernel : public ApplyMonteCarloBarostatKernel {
public:
1021
    ReferenceApplyMonteCarloBarostatKernel(std::string name, const Platform& platform) : ApplyMonteCarloBarostatKernel(name, platform), barostat(NULL) {
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
    }
    ~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;
};

1052
1053
1054
1055
1056
/**
 * This kernel is invoked to calculate the kinetic energy of the system.
 */
class ReferenceCalcKineticEnergyKernel : public CalcKineticEnergyKernel {
public:
1057
    ReferenceCalcKineticEnergyKernel(std::string name, const Platform& platform) : CalcKineticEnergyKernel(name, platform) {
1058
1059
    }
    /**
1060
     * Initialize the kernel.
1061
     * 
1062
     * @param system     the System this kernel will be applied to
1063
     */
1064
    void initialize(const System& system);
1065
1066
1067
    /**
     * Execute the kernel.
     * 
1068
     * @param context    the context in which to execute this kernel
1069
     */
1070
    double execute(ContextImpl& context);
1071
1072
private:
    std::vector<double> masses;
1073
1074
};

1075
1076
1077
1078
1079
/**
 * This kernel is invoked to remove center of mass motion from the system.
 */
class ReferenceRemoveCMMotionKernel : public RemoveCMMotionKernel {
public:
1080
    ReferenceRemoveCMMotionKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : RemoveCMMotionKernel(name, platform), data(data) {
1081
1082
    }
    /**
Peter Eastman's avatar
Peter Eastman committed
1083
     * Initialize the kernel, setting up the particle masses.
1084
     * 
1085
1086
     * @param system     the System this kernel will be applied to
     * @param force      the CMMotionRemover this kernel will be used for
1087
     */
1088
    void initialize(const System& system, const CMMotionRemover& force);
1089
1090
1091
    /**
     * Execute the kernel.
     * 
1092
     * @param context    the context in which to execute this kernel
1093
     */
1094
    void execute(ContextImpl& context);
1095
private:
1096
    ReferencePlatform::PlatformData& data;
1097
    std::vector<double> masses;
1098
    int frequency;
1099
1100
};

1101
1102
1103
} // namespace OpenMM

#endif /*OPENMM_REFERENCEKERNELS_H_*/