ReferenceKernels.h 84.2 KB
Newer Older
1
2
3
4
5
6
#ifndef OPENMM_REFERENCEKERNELS_H_
#define OPENMM_REFERENCEKERNELS_H_

/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
Evan Pretti's avatar
Evan Pretti committed
7
8
 * This is part of the OpenMM molecular simulation toolkit.                   *
 * See https://openmm.org/development.                                        *
9
 *                                                                            *
10
 * Portions copyright (c) 2008-2025 Stanford University and the Authors.      *
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 * 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.                                     *
 * -------------------------------------------------------------------------- */

33
#include "ReferencePlatform.h"
34
#include "openmm/kernels.h"
35
#include "openmm/internal/CustomCPPForceImpl.h"
36
#include "openmm/internal/CustomNonbondedForceImpl.h"
37
#include "openmm/internal/windowsExport.h"
38
#include "SimTKOpenMMRealType.h"
39
#include "ReferenceConstantPotential.h"
40
#include "ReferenceNeighborList.h"
41
#include "lepton/CompiledExpression.h"
42
#include "lepton/CustomFunction.h"
43
44
#include <array>
#include <utility>
45

46
47
namespace OpenMM {

48
class ReferenceObc;
49
class ReferenceAndersenThermostat;
50
class ReferenceLangevinMiddleDynamics;
51
52
53
54
class ReferenceCustomBondIxn;
class ReferenceCustomAngleIxn;
class ReferenceCustomTorsionIxn;
class ReferenceCustomExternalIxn;
55
class ReferenceCustomCentroidBondIxn;
56
class ReferenceCustomCompoundBondIxn;
57
class ReferenceCustomCVForce;
58
class ReferenceCustomHbondIxn;
59
class ReferenceCustomManyParticleIxn;
60
class ReferenceGayBerneForce;
61
class ReferenceBrownianDynamics;
62
class ReferenceConstraintAlgorithm;
63
class ReferenceNoseHooverChain;
64
class ReferenceMonteCarloBarostat;
65
class ReferenceNoseHooverDynamics;
66
class ReferenceVariableStochasticDynamics;
67
class ReferenceVariableVerletDynamics;
68
class ReferenceVerletDynamics;
69
class ReferenceCustomDynamics;
Peter Eastman's avatar
Peter Eastman committed
70
class ReferenceDPDDynamics;
71
class ReferenceQTBDynamics;
72

73
/**
74
75
76
 * 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.
77
 */
78
class ReferenceCalcForcesAndEnergyKernel : public CalcForcesAndEnergyKernel {
79
public:
80
    ReferenceCalcForcesAndEnergyKernel(std::string name, const Platform& platform) : CalcForcesAndEnergyKernel(name, platform) {
81
82
83
84
85
86
87
88
    }
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     */
    void initialize(const System& system);
    /**
89
     * This is called at the beginning of each force/energy computation, before calcForcesAndEnergy() has been called on
90
91
     * any ForceImpl.
     *
92
93
94
     * @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
95
     * @param groups        a set of bit flags for which force groups to include
96
     */
97
    void beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups);
98
    /**
99
     * This is called at the end of each force/energy computation, after calcForcesAndEnergy() has been called on
100
101
     * every ForceImpl.
     *
102
103
104
     * @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
105
     * @param groups        a set of bit flags for which force groups to include
106
107
     * @param valid         the method may set this to false to indicate the results are invalid and the force/energy
     *                      calculation should be repeated
108
     * @return the potential energy of the system.  This value is added to all values returned by ForceImpls'
109
     * calcForcesAndEnergy() methods.  That is, each force kernel may <i>either</i> return its contribution to the
110
     * energy directly, <i>or</i> add it to an internal buffer so that it will be included here.
111
     */
112
    double finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups, bool& valid);
113
private:
peastman's avatar
peastman committed
114
    std::vector<Vec3> savedForces;
115
116
};

117
/**
118
119
 * This kernel provides methods for setting and retrieving various state data: time, positions,
 * velocities, and forces.
120
 */
121
class OPENMM_EXPORT ReferenceUpdateStateDataKernel : public UpdateStateDataKernel {
122
public:
123
    ReferenceUpdateStateDataKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : UpdateStateDataKernel(name, platform), data(data) {
124
125
126
127
128
129
130
131
132
133
134
135
    }
    /**
     * 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
     */
136
    double getTime(const ContextImpl& context) const;
137
138
139
140
141
    /**
     * Set the current time (in picoseconds).
     *
     * @param context    the context in which to execute this kernel
     */
142
    void setTime(ContextImpl& context, double time);
143
144
145
146
147
148
149
150
151
152
153
154
    /**
     * Get the current step count
     *
     * @param context    the context in which to execute this kernel
     */
    long long getStepCount(const ContextImpl& context) const;
    /**
     * Set the current step count
     *
     * @param context    the context in which to execute this kernel
     */
    void setStepCount(const ContextImpl& context, long long count);
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
    /**
     * 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);
179
180
181
182
183
184
185
186
187
    /**
     * Compute velocities, shifted in time to account for a leapfrog integrator.  The shift
     * is based on the most recently computed forces.
     * 
     * @param context     the context in which to execute this kernel
     * @param timeShift   the amount by which to shift the velocities in time
     * @param velocities  the shifted velocities are returned in this
     */
    void computeShiftedVelocities(ContextImpl& context, double timeShift, std::vector<Vec3>& velocities);
188
189
190
191
192
193
    /**
     * Get the current forces on all particles.
     *
     * @param forces  on exit, this contains the forces
     */
    void getForces(ContextImpl& context, std::vector<Vec3>& forces);
194
195
196
197
198
199
    /**
     * Get the current derivatives of the energy with respect to context parameters.
     *
     * @param derivs  on exit, this contains the derivatives
     */
    void getEnergyParameterDerivatives(ContextImpl& context, std::map<std::string, double>& derivs);
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
    /**
     * 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
     */
215
    void setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c);
Peter Eastman's avatar
Peter Eastman committed
216
217
218
219
220
221
222
223
224
225
226
227
    /**
     * 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);
228
229
private:
    ReferencePlatform::PlatformData& data;
230
    std::vector<double> masses;
231
232
};

233
234
235
236
237
238
/**
 * 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) :
239
            ApplyConstraintsKernel(name, platform), data(data) {
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
    }
    ~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);
255
256
257
258
259
260
261
    /**
     * Update particle velocities to enforce constraints.
     *
     * @param context    the context in which to execute this kernel
     * @param tol        the velocity tolerance within which constraints must be satisfied.
     */
    void applyToVelocities(ContextImpl& context, double tol);
262
263
private:
    ReferencePlatform::PlatformData& data;
peastman's avatar
peastman committed
264
265
    std::vector<double> masses;
    std::vector<double> inverseMasses;
266
267
};

268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
/**
 * 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);
};

289
/**
290
 * This kernel is invoked by HarmonicBondForce to calculate the forces acting on the system and the energy of the system.
291
 */
292
class ReferenceCalcHarmonicBondForceKernel : public CalcHarmonicBondForceKernel {
293
public:
294
    ReferenceCalcHarmonicBondForceKernel(std::string name, const Platform& platform) : CalcHarmonicBondForceKernel(name, platform) {
295
296
    }
    /**
297
     * Initialize the kernel.
298
     * 
299
     * @param system     the System this kernel will be applied to
300
301
302
303
     * @param force      the HarmonicBondForce this kernel will be used for
     */
    void initialize(const System& system, const HarmonicBondForce& force);
    /**
304
305
306
307
308
309
     * 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
310
     */
311
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
312
313
314
315
316
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the HarmonicBondForce to copy the parameters from
317
318
     * @param firstBond  the index of the first bond whose parameters might have changed
     * @param lastBond   the index of the last bond whose parameters might have changed
319
     */
320
    void copyParametersToContext(ContextImpl& context, const HarmonicBondForce& force, int firstBond, int lastBond);
321
322
private:
    int numBonds;
323
324
    std::vector<std::vector<int> >bondIndexArray;
    std::vector<std::vector<double> >bondParamArray;
325
    bool usePeriodic;
326
327
};

328
329
330
331
332
/**
 * 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:
333
    ReferenceCalcCustomBondForceKernel(std::string name, const Platform& platform) : CalcCustomBondForceKernel(name, platform), ixn(NULL) {
334
    }
335
    ~ReferenceCalcCustomBondForceKernel();
336
337
338
339
340
341
342
343
    /**
     * 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);
    /**
344
     * Execute the kernel to calculate the forces and/or energy.
345
     *
346
347
348
349
     * @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
350
     */
351
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
352
353
354
    /**
     * Copy changed parameters over to a context.
     *
355
356
     * @param firstBond  the index of the first bond whose parameters might have changed
     * @param lastBond   the index of the last bond whose parameters might have changed
357
     */
358
    void copyParametersToContext(ContextImpl& context, const CustomBondForce& force, int firstBond, int lastBond);
359
360
private:
    int numBonds;
361
    ReferenceCustomBondIxn* ixn;
362
363
    std::vector<std::vector<int> >bondIndexArray;
    std::vector<std::vector<double> >bondParamArray;
364
    Lepton::CompiledExpression energyExpression, forceExpression;
365
366
    std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
    std::vector<std::string> parameterNames, globalParameterNames, energyParamDerivNames;
367
    bool usePeriodic;
368
369
};

370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
/**
 * 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) {
    }
    /**
     * 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);
    /**
385
386
387
388
389
390
     * 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
391
     */
392
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
393
394
395
    /**
     * Copy changed parameters over to a context.
     *
396
397
     * @param firstAngle the index of the first bond whose parameters might have changed
     * @param lastAngle  the index of the last bond whose parameters might have changed
398
     */
399
    void copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force, int firstAngle, int lastAngle);
400
401
private:
    int numAngles;
402
403
    std::vector<std::vector<int> >angleIndexArray;
    std::vector<std::vector<double> >angleParamArray;
404
    bool usePeriodic;
405
406
};

407
408
409
410
411
/**
 * 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:
412
    ReferenceCalcCustomAngleForceKernel(std::string name, const Platform& platform) : CalcCustomAngleForceKernel(name, platform), ixn(NULL) {
413
    }
414
    ~ReferenceCalcCustomAngleForceKernel();
415
416
417
418
419
420
421
422
    /**
     * 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);
    /**
423
     * Execute the kernel to calculate the forces and/or energy.
424
     *
425
426
427
428
     * @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
429
     */
430
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
431
432
433
    /**
     * Copy changed parameters over to a context.
     *
434
435
     * @param firstAngle the index of the first bond whose parameters might have changed
     * @param lastAngle  the index of the last bond whose parameters might have changed
436
     */
437
    void copyParametersToContext(ContextImpl& context, const CustomAngleForce& force, int firstAngle, int lastAngle);
438
439
private:
    int numAngles;
440
    ReferenceCustomAngleIxn* ixn;
441
442
    std::vector<std::vector<int> >angleIndexArray;
    std::vector<std::vector<double> >angleParamArray;
443
    Lepton::CompiledExpression energyExpression, forceExpression;
444
445
    std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
    std::vector<std::string> parameterNames, globalParameterNames, energyParamDerivNames;
446
    bool usePeriodic;
447
448
};

449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
/**
 * 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) {
    }
    /**
     * 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);
    /**
464
465
466
467
468
469
     * 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
470
     */
471
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
472
473
474
    /**
     * Copy changed parameters over to a context.
     *
475
476
477
478
     * @param context      the context to copy parameters to
     * @param force        the PeriodicTorsionForce to copy the parameters from
     * @param firstTorsion the index of the first torsion whose parameters might have changed
     * @param lastTorsion  the index of the last torsion whose parameters might have changed
479
     */
480
    void copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force, int firstTorsion, int lastTorsion);
481
482
private:
    int numTorsions;
483
484
    std::vector<std::vector<int> >torsionIndexArray;
    std::vector<std::vector<double> >torsionParamArray;
485
    bool usePeriodic;
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
};

/**
 * 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) {
    }
    /**
     * 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);
    /**
503
504
505
506
507
508
     * 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
509
     */
510
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
511
512
513
514
515
516
517
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the RBTorsionForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const RBTorsionForce& force);
518
519
private:
    int numTorsions;
520
521
    std::vector<std::vector<int> >torsionIndexArray;
    std::vector<std::vector<double> >torsionParamArray;
522
    bool usePeriodic;
523
524
};

525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
/**
 * 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);
    /**
540
     * Execute the kernel to calculate the forces and/or energy.
541
     *
542
543
544
545
     * @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
546
     */
547
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
548
549
550
551
552
553
554
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CMAPTorsionForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force);
555
private:
peastman's avatar
peastman committed
556
    std::vector<std::vector<std::vector<double> > > coeff;
557
558
    std::vector<int> torsionMaps;
    std::vector<std::vector<int> > torsionIndices;
559
    bool usePeriodic;
560
561
};

562
563
564
565
566
/**
 * 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:
567
    ReferenceCalcCustomTorsionForceKernel(std::string name, const Platform& platform) : CalcCustomTorsionForceKernel(name, platform), ixn(NULL) {
568
    }
569
    ~ReferenceCalcCustomTorsionForceKernel();
570
571
572
573
574
575
576
577
    /**
     * 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);
    /**
578
     * Execute the kernel to calculate the forces and/or energy.
579
     *
580
581
582
583
     * @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
584
     */
585
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
586
587
588
    /**
     * Copy changed parameters over to a context.
     *
589
590
591
592
     * @param context      the context to copy parameters to
     * @param force        the CustomTorsionForce to copy the parameters from
     * @param firstTorsion the index of the first torsion whose parameters might have changed
     * @param lastTorsion  the index of the last torsion whose parameters might have changed
593
     */
594
    void copyParametersToContext(ContextImpl& context, const CustomTorsionForce& force, int firstTorsion, int lastTorsion);
595
596
private:
    int numTorsions;
597
    ReferenceCustomTorsionIxn* ixn;
598
599
    std::vector<std::vector<int> >torsionIndexArray;
    std::vector<std::vector<double> >torsionParamArray;
600
    Lepton::CompiledExpression energyExpression, forceExpression;
601
602
    std::vector<Lepton::CompiledExpression> energyParamDerivExpressions;
    std::vector<std::string> parameterNames, globalParameterNames, energyParamDerivNames;
603
    bool usePeriodic;
604
605
};

606
607
608
609
610
611
612
613
614
615
616
617
618
/**
 * 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
619
     */
620
    void initialize(const System& system, const NonbondedForce& force);
621
    /**
622
623
624
625
626
     * 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
627
     * @param includeReciprocal  true if reciprocal space interactions should be included
628
     * @return the potential energy due to the force
629
     */
630
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal);
631
632
633
    /**
     * Copy changed parameters over to a context.
     *
634
635
636
637
638
639
     * @param context        the context to copy parameters to
     * @param force          the NonbondedForce to copy the parameters from
     * @param firstParticle  the index of the first particle whose parameters might have changed
     * @param lastParticle   the index of the last particle whose parameters might have changed
     * @param firstException the index of the first exception whose parameters might have changed
     * @param lastException  the index of the last exception whose parameters might have changed
640
     */
641
    void copyParametersToContext(ContextImpl& context, const NonbondedForce& force, int firstParticle, int lastParticle, int firstException, int lastException);
642
643
644
645
646
647
648
649
650
    /**
     * Get the parameters being used for PME.
     * 
     * @param alpha   the separation parameter
     * @param nx      the number of grid points along the X axis
     * @param ny      the number of grid points along the Y axis
     * @param nz      the number of grid points along the Z axis
     */
    void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
651
652
653
    /**
     * Get the dispersion parameters being used for the dispersion term in LJPME.
     *
654
655
656
657
     * @param alpha   the separation parameter
     * @param nx      the number of grid points along the X axis
     * @param ny      the number of grid points along the Y axis
     * @param nz      the number of grid points along the Z axis
658
     */
659
    void getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
660
private:
661
    void computeParameters(ContextImpl& context);
Peter Eastman's avatar
Peter Eastman committed
662
    int numParticles, num14;
663
664
    std::vector<std::vector<int> >bonded14IndexArray;
    std::vector<std::vector<double> > particleParamArray, bonded14ParamArray;
665
666
    std::vector<std::array<double, 3> > baseParticleParams, baseExceptionParams;
    std::map<std::pair<std::string, int>, std::array<double, 3> > particleParamOffsets, exceptionParamOffsets;
667
    std::map<int, int> nb14Index;
peastman's avatar
peastman committed
668
    double nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, ewaldDispersionAlpha, dispersionCoefficient;
669
    int kmax[3], gridSize[3], dispersionGridSize[3];
670
    bool useSwitchingFunction, exceptionsArePeriodic;
671
672
673
    std::vector<std::set<int> > exclusions;
    NonbondedMethod nonbondedMethod;
    NeighborList* neighborList;
674
675
};

676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
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
747
/**
 * This kernel is invoked by ConstantPotentialForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcConstantPotentialForceKernel : public CalcConstantPotentialForceKernel {
public:
    ReferenceCalcConstantPotentialForceKernel(std::string name, const Platform& platform) : CalcConstantPotentialForceKernel(name, platform), neighborList(NULL), solver(NULL) {
    }
    ~ReferenceCalcConstantPotentialForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param force      the ConstantPotentialForce this kernel will be used for
     */
    void initialize(const System& system, const ConstantPotentialForce& force);
    /**
     * Execute the kernel to calculate the forces and/or energy.
     *
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
     */
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
    /**
     * Copy changed parameters over to a context.
     *
     * @param context        the context to copy parameters to
     * @param force          the ConstantPotentialForce to copy the parameters from
     * @param firstParticle  the index of the first particle whose parameters might have changed
     * @param lastParticle   the index of the last particle whose parameters might have changed
     * @param firstException the index of the first exception whose parameters might have changed
     * @param lastException  the index of the last exception whose parameters might have changed
     * @param firstElectrode the index of the first electrode whose parameters might have changed
     * @param lastElectrode  the index of the last electrode whose parameters might have changed
     */
    void copyParametersToContext(ContextImpl& context, const ConstantPotentialForce& force, int firstParticle, int lastParticle, int firstException, int lastException, int firstElectrode, int lastElectrode);
    /**
     * Get the parameters being used for PME.
     *
     * @param alpha   the separation parameter
     * @param nx      the number of grid points along the X axis
     * @param ny      the number of grid points along the Y axis
     * @param nz      the number of grid points along the Z axis
     */
    void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
    /**
     * Get the charges on all particles.
     *
     * @param context       the context to copy parameters to
     * @param[out] charges  a vector to populate with particle charges
     */
    void getCharges(ContextImpl& context, std::vector<double>& charges);
private:
    void updateNeighborList(const Vec3* boxVectors, const std::vector<Vec3>& posData);
private:
    int numParticles, num14, numElectrodeParticles;
    std::vector<double> charges;
    std::vector<std::vector<double> > bonded14ParamArray;
    std::vector<std::vector<int> > bonded14IndexArray;
    std::map<int, int> nb14Index;
    std::vector<std::set<int> > exclusions;
    std::vector<int> sysToElec, elecToSys;
    std::vector<std::array<double, 3> > electrodeParamArray;
    double nonbondedCutoff, ewaldAlpha, cgErrorTol, chargeTarget;
    int gridSize[3];
    bool exceptionsArePeriodic, useChargeConstraint;
    Vec3 externalField;
    NeighborList* neighborList;
    ReferenceConstantPotentialSolver* solver;
};

748
749
750
751
752
/**
 * This kernel is invoked by CustomNonbondedForce to calculate the forces acting on the system.
 */
class ReferenceCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
public:
753
    ReferenceCalcCustomNonbondedForceKernel(std::string name, const Platform& platform) : CalcCustomNonbondedForceKernel(name, platform), forceCopy(NULL) {
754
755
756
757
758
759
760
761
762
763
    }
    ~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);
    /**
764
     * Execute the kernel to calculate the forces and/or energy.
765
     *
766
767
768
769
     * @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
770
     */
771
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
772
773
774
775
776
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomNonbondedForce to copy the parameters from
777
778
     * @param firstParticle  the index of the first particle whose parameters might have changed
     * @param lastParticle   the index of the last particle whose parameters might have changed
779
     */
780
    void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force, int firstParticle, int lastParticle);
781
private:
782
    void createExpressions(const CustomNonbondedForce& force);
783
    int numParticles;
784
    std::vector<std::vector<double> > particleParamArray;
peastman's avatar
peastman committed
785
    double nonbondedCutoff, switchingDistance, periodicBoxSize[3], longRangeCoefficient;
786
    bool useSwitchingFunction, hasInitializedLongRangeCorrection;
787
    CustomNonbondedForce* forceCopy;
788
    CustomNonbondedForceImpl::LongRangeCorrectionData longRangeCorrectionData;
789
    std::map<std::string, double> globalParamValues;
790
    std::vector<std::set<int> > exclusions;
791
    Lepton::CompiledExpression energyExpression, forceExpression;
792
793
    std::vector<Lepton::CompiledExpression> computedValueExpressions, energyParamDerivExpressions;
    std::vector<std::string> parameterNames, globalParameterNames, computedValueNames, energyParamDerivNames;
794
    std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
795
    std::vector<double> longRangeCoefficientDerivs;
796
    std::map<std::string, int> tabulatedFunctionUpdateCount;
797
798
799
800
    NonbondedMethod nonbondedMethod;
    NeighborList* neighborList;
};

801
/**
802
 * This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
803
 */
804
class ReferenceCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
805
public:
806
    ReferenceCalcGBSAOBCForceKernel(std::string name, const Platform& platform) : CalcGBSAOBCForceKernel(name, platform) {
807
    }
808
    ~ReferenceCalcGBSAOBCForceKernel();
809
    /**
810
     * Initialize the kernel.
811
     * 
812
     * @param system     the System this kernel will be applied to
813
     * @param force      the GBSAOBCForce this kernel will be used for
814
     */
815
    void initialize(const System& system, const GBSAOBCForce& force);
816
    /**
817
818
819
820
821
822
     * 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
823
     */
824
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
825
826
827
828
829
830
831
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the GBSAOBCForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
832
private:
833
    ReferenceObc* obc;
peastman's avatar
peastman committed
834
    std::vector<double> charges;
835
    bool isPeriodic;
836
837
};

838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
/**
 * 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);
    /**
854
     * Execute the kernel to calculate the forces and/or energy.
855
     *
856
857
858
859
     * @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
860
     */
861
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
862
863
864
865
866
867
868
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomGBForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomGBForce& force);
869
private:
870
    void createExpressions(const CustomGBForce& force);
871
    int numParticles;
872
    bool isPeriodic;
873
    std::vector<std::vector<double> > particleParamArray;
peastman's avatar
peastman committed
874
    double nonbondedCutoff;
875
    std::vector<std::set<int> > exclusions;
876
    std::vector<std::string> particleParameterNames, globalParameterNames, energyParamDerivNames, valueNames;
877
878
879
    std::vector<Lepton::CompiledExpression> valueExpressions;
    std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions;
    std::vector<std::vector<Lepton::CompiledExpression> > valueGradientExpressions;
880
    std::vector<std::vector<Lepton::CompiledExpression> > valueParamDerivExpressions;
881
    std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
882
883
884
    std::vector<Lepton::CompiledExpression> energyExpressions;
    std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions;
    std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions;
885
    std::vector<std::vector<Lepton::CompiledExpression> > energyParamDerivExpressions;
886
    std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
887
    std::map<std::string, int> tabulatedFunctionUpdateCount;
888
889
890
891
    NonbondedMethod nonbondedMethod;
    NeighborList* neighborList;
};

892
893
894
895
896
/**
 * 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:
897
    ReferenceCalcCustomExternalForceKernel(std::string name, const Platform& platform) : CalcCustomExternalForceKernel(name, platform), ixn(NULL) {
898
    }
899
    ~ReferenceCalcCustomExternalForceKernel();
900
901
902
903
904
905
906
907
    /**
     * 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);
    /**
908
     * Execute the kernel to calculate the forces and/or energy.
909
     *
910
911
912
913
     * @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
914
     */
915
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
916
917
918
    /**
     * Copy changed parameters over to a context.
     *
919
920
921
922
     * @param context        the context to copy parameters to
     * @param force          the CustomExternalForce to copy the parameters from
     * @param firstParticle  the index of the first particle whose parameters might have changed
     * @param lastParticle   the index of the last particle whose parameters might have changed
923
     */
924
    void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force, int firstParticle, int lastParticle);
925
926
private:
    int numParticles;
927
    ReferenceCustomExternalIxn* ixn;
928
    std::vector<int> particles;
929
    std::vector<std::vector<double> > particleParamArray;
930
    Lepton::CompiledExpression energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ;
931
    std::vector<std::string> parameterNames, globalParameterNames;
peastman's avatar
peastman committed
932
    Vec3* boxVectors;
933
934
};

935
936
937
938
939
/**
 * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system.
 */
class ReferenceCalcCustomHbondForceKernel : public CalcCustomHbondForceKernel {
public:
940
    ReferenceCalcCustomHbondForceKernel(std::string name, const Platform& platform) : CalcCustomHbondForceKernel(name, platform), ixn(NULL) {
941
942
943
944
945
946
947
948
949
950
    }
    ~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);
    /**
951
     * Execute the kernel to calculate the forces and/or energy.
952
     *
953
954
955
956
     * @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
957
     */
958
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
959
960
961
962
963
964
965
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomHbondForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomHbondForce& force);
966
private:
967
    void createInteraction(const CustomHbondForce& force);
968
    int numDonors, numAcceptors, numParticles;
969
    bool isPeriodic;
970
    std::vector<std::vector<int> > donorParticles, acceptorParticles;
971
    std::vector<std::vector<double> > donorParamArray, acceptorParamArray;
peastman's avatar
peastman committed
972
    double nonbondedCutoff;
973
    ReferenceCustomHbondIxn* ixn;
974
    std::vector<std::set<int> > exclusions;
975
    std::vector<std::string> globalParameterNames;
976
    std::map<std::string, int> tabulatedFunctionUpdateCount;
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
1008
1009
1010
/**
 * This kernel is invoked by CustomCentroidBondForce to calculate the forces acting on the system.
 */
class ReferenceCalcCustomCentroidBondForceKernel : public CalcCustomCentroidBondForceKernel {
public:
    ReferenceCalcCustomCentroidBondForceKernel(std::string name, const Platform& platform) : CalcCustomCentroidBondForceKernel(name, platform), ixn(NULL) {
    }
    ~ReferenceCalcCustomCentroidBondForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomCentroidBondForce this kernel will be used for
     */
    void initialize(const System& system, const CustomCentroidBondForce& force);
    /**
     * Execute the kernel to calculate the forces and/or energy.
     *
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
     */
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomCentroidBondForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force);
private:
1011
    void createInteraction(const CustomCentroidBondForce& force);
1012
    int numBonds, numParticles;
1013
1014
1015
    std::vector<std::vector<int> > bondGroups;
    std::vector<std::vector<int> > groupAtoms;
    std::vector<std::vector<double> > normalizedWeights;
1016
    std::vector<std::vector<double> > bondParamArray;
1017
    ReferenceCustomCentroidBondIxn* ixn;
1018
    std::vector<std::string> globalParameterNames, energyParamDerivNames;
1019
    std::map<std::string, int> tabulatedFunctionUpdateCount;
1020
    bool usePeriodic;
1021
    Vec3* boxVectors;
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
/**
 * 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);
1048
1049
1050
1051
1052
1053
1054
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomCompoundBondForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force);
1055
private:
1056
    void createInteraction(const CustomCompoundBondForce& force);
1057
    int numBonds;
1058
    std::vector<std::vector<int> > bondParticles;
1059
    std::vector<std::vector<double> > bondParamArray;
1060
    ReferenceCustomCompoundBondIxn* ixn;
1061
    std::vector<std::string> globalParameterNames, energyParamDerivNames;
1062
    std::map<std::string, int> tabulatedFunctionUpdateCount;
1063
    bool usePeriodic;
1064
    Vec3* boxVectors;
1065
1066
};

1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
/**
 * This kernel is invoked by CustomManyParticleForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcCustomManyParticleForceKernel : public CalcCustomManyParticleForceKernel {
public:
    ReferenceCalcCustomManyParticleForceKernel(std::string name, const Platform& platform) : CalcCustomManyParticleForceKernel(name, platform), ixn(NULL) {
    }
    ~ReferenceCalcCustomManyParticleForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomManyParticleForce this kernel will be used for
     */
    void initialize(const System& system, const CustomManyParticleForce& force);
    /**
     * Execute the kernel to calculate the forces and/or energy.
     *
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
     */
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomManyParticleForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force);
private:
    int numParticles;
peastman's avatar
peastman committed
1100
    double cutoffDistance;
1101
    std::vector<std::vector<double> > particleParamArray;
1102
1103
    ReferenceCustomManyParticleIxn* ixn;
    std::vector<std::string> globalParameterNames;
1104
    std::map<std::string, int> tabulatedFunctionUpdateCount;
1105
1106
1107
    NonbondedMethod nonbondedMethod;
};

1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
/**
 * This kernel is invoked by GayBerneForce to calculate the forces acting on the system.
 */
class ReferenceCalcGayBerneForceKernel : public CalcGayBerneForceKernel {
public:
    ReferenceCalcGayBerneForceKernel(std::string name, const Platform& platform) : CalcGayBerneForceKernel(name, platform), ixn(NULL) {
    }
    ~ReferenceCalcGayBerneForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the GayBerneForce this kernel will be used for
     */
    void initialize(const System& system, const GayBerneForce& force);
    /**
     * Execute the kernel to calculate the forces and/or energy.
     *
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @return the potential energy due to the force
     */
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the GayBerneForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const GayBerneForce& force);
private:
    ReferenceGayBerneForce* ixn;
};

1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
/**
 * This kernel is invoked by CustomCVForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcCustomCVForceKernel : public CalcCustomCVForceKernel {
public:
    ReferenceCalcCustomCVForceKernel(std::string name, const Platform& platform) : CalcCustomCVForceKernel(name, platform), ixn(NULL) {
    }
    ~ReferenceCalcCustomCVForceKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomCVForce this kernel will be used for
1155
     * @param innerContext   the context created by the CustomCVForce for computing collective variables
1156
     */
1157
    void initialize(const System& system, const CustomCVForce& force, ContextImpl& innerContext);
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
    /**
     * Execute the kernel to calculate the forces and/or energy.
     *
     * @param context        the context in which to execute this kernel
     * @param innerContext   the context created by the CustomCVForce for computing collective variables
     * @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, ContextImpl& innerContext, bool includeForces, bool includeEnergy);
    /**
     * Copy state information to the inner context.
     *
     * @param context        the context in which to execute this kernel
     * @param innerContext   the context created by the CustomCVForce for computing collective variables
     */
    void copyState(ContextImpl& context, ContextImpl& innerContext);
1175
1176
1177
1178
1179
1180
1181
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the CustomCVForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const CustomCVForce& force);
1182
1183
1184
1185
1186
private:
    ReferenceCustomCVForce* ixn;
    std::vector<std::string> globalParameterNames, energyParamDerivNames;
};

1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
/**
 * This kernel is invoked by RMSDForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcRMSDForceKernel : public CalcRMSDForceKernel {
public:
    ReferenceCalcRMSDForceKernel(std::string name, const Platform& platform) : CalcRMSDForceKernel(name, platform) {
    }
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the RMSDForce this kernel will be used for
     */
    void initialize(const System& system, const RMSDForce& force);
    /**
     * Execute the kernel to calculate the forces and/or energy.
     *
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
     */
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the RMSDForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const RMSDForce& force);
private:
    std::vector<Vec3> referencePos;
    std::vector<int> particles;
};

1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
/**
 * This kernel is invoked by RGForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcRGForceKernel : public CalcRGForceKernel {
public:
    ReferenceCalcRGForceKernel(std::string name, const Platform& platform) : CalcRGForceKernel(name, platform) {
    }
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the RGForce this kernel will be used for
     */
    void initialize(const System& system, const RGForce& 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:
    std::vector<int> particles;
};

1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
/**
 * This kernel is invoked by OrientationRestraintForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcOrientationRestraintForceKernel : public CalcOrientationRestraintForceKernel {
public:
    ReferenceCalcOrientationRestraintForceKernel(std::string name, const Platform& platform) : CalcOrientationRestraintForceKernel(name, platform) {
    }
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the OrientationRestraintForce this kernel will be used for
     */
    void initialize(const System& system, const OrientationRestraintForce& force);
    /**
     * Execute the kernel to calculate the forces and/or energy.
     *
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
     */
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the OrientationRestraintForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const OrientationRestraintForce& force);
private:
    double k;
    std::vector<Vec3> referencePos;
    std::vector<int> particles;
};

1285
1286
1287
1288
1289
/**
 * This kernel is invoked by VerletIntegrator to take one time step.
 */
class ReferenceIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
public:
1290
    ReferenceIntegrateVerletStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVerletStepKernel(name, platform),
1291
        data(data), dynamics(0) {
1292
    }
1293
    ~ReferenceIntegrateVerletStepKernel();
1294
    /**
1295
     * Initialize the kernel.
1296
     * 
1297
1298
     * @param system     the System this kernel will be applied to
     * @param integrator the VerletIntegrator this kernel will be used for
1299
     */
1300
    void initialize(const System& system, const VerletIntegrator& integrator);
1301
1302
1303
    /**
     * Execute the kernel.
     * 
1304
1305
     * @param context    the context in which to execute this kernel
     * @param integrator the VerletIntegrator this kernel is being used for
1306
     */
1307
    void execute(ContextImpl& context, const VerletIntegrator& integrator);
1308
1309
1310
1311
1312
1313
1314
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the VerletIntegrator this kernel is being used for
     */
    double computeKineticEnergy(ContextImpl& context, const VerletIntegrator& integrator);
1315
private:
1316
    ReferencePlatform::PlatformData& data;
1317
    ReferenceVerletDynamics* dynamics;
peastman's avatar
peastman committed
1318
    std::vector<double> masses;
1319
    double prevStepSize;
1320
1321
1322
}
;
/**
1323
 * This kernel is invoked by NoseHooverIntegrator to take one time step.
1324
 */
1325
class ReferenceIntegrateNoseHooverStepKernel : public IntegrateNoseHooverStepKernel {
1326
public:
1327
    ReferenceIntegrateNoseHooverStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateNoseHooverStepKernel(name, platform),
1328
1329
        data(data), dynamics(0) {
    }
1330
    ~ReferenceIntegrateNoseHooverStepKernel();
1331
1332
1333
1334
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
1335
     * @param integrator the NoseHooverIntegrator this kernel will be used for
1336
     */
1337
    void initialize(const System& system, const NoseHooverIntegrator& integrator);
1338
1339
1340
1341
1342
1343
    /**
     * Execute the kernel.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the VerletIntegrator this kernel is being used for
     */
1344
    void execute(ContextImpl& context, const NoseHooverIntegrator& integrator);
1345
1346
1347
1348
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
1349
     * @param integrator the NoseHooverIntegrator this kernel is being used for
1350
     */
1351
    double computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator);
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
    /**
     * Execute the kernel that propagates the Nose Hoover chain and determines the velocity scale factor.
     * 
     * @param context  the context in which to execute this kernel
     * @param noseHooverChain the object describing the chain to be propagated.
     * @param kineticEnergy the {center of mass, relative} kineticEnergies of the particles being thermostated by this chain.
     * @param timeStep the time step used by the integrator.
     * @return the velocity scale factor to apply to the particles associated with this heat bath.
     */
    std::pair<double, double> propagateChain(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> kineticEnergy, double timeStep);
    /**
     * Execute the kernal that computes the total (kinetic + potential) heat bath energy.
     *
     * @param context the context in which to execute this kernel
     * @param noseHooverChain the chain whose energy is to be determined.
     * @return the total heat bath energy.
     */
    double computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain);
    /**
     * Execute the kernel that computes the kinetic energy for a subset of atoms,
     * or the relative kinetic energy of Drude particles with respect to their parent atoms
     *
     * @param context the context in which to execute this kernel
     * @param noseHooverChain the chain whose energy is to be determined.
     * @param downloadValue whether the computed value should be downloaded and returned.
     */
    std::pair<double, double> computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue);
    /**
     * Execute the kernel that scales the velocities of particles associated with a nose hoover chain
     *
     * @param context the context in which to execute this kernel
     * @param noseHooverChain the chain whose energy is to be determined.
     * @param scaleFactor the multiplicative factor by which {absolute, relative} velocities are scaled.
     */
    void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactor);
1387
1388
1389
1390
1391
1392
1393
1394
    /**
     * Write the chain states to a checkpoint.
     */
    void createCheckpoint(ContextImpl& context, std::ostream& stream) const;
    /**
     * Load the chain states from a checkpoint.
     */
    void loadCheckpoint(ContextImpl& context, std::istream& stream);
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
    /**
     * Get the internal states of all chains.
     * 
     * @param context       the context for which to get the states
     * @param positions     element [i][j] contains the position of bead j for chain i
     * @param velocities    element [i][j] contains the velocity of bead j for chain i
     */
    void getChainStates(ContextImpl& context, std::vector<std::vector<double> >& positions, std::vector<std::vector<double> >& velocities) const;
    /**
     * Set the internal states of all chains.
     * 
     * @param context       the context for which to get the states
     * @param positions     element [i][j] contains the position of bead j for chain i
     * @param velocities    element [i][j] contains the velocity of bead j for chain i
     */
    void setChainStates(ContextImpl& context, const std::vector<std::vector<double> >& positions, const std::vector<std::vector<double> >& velocities);
1411
1412
private:
    ReferencePlatform::PlatformData& data;
1413
1414
    ReferenceNoseHooverChain* chainPropagator;
    ReferenceNoseHooverDynamics* dynamics;
1415
    std::vector<double> masses;
1416
1417
    std::vector<std::vector<double> > chainPositions;
    std::vector<std::vector<double> > chainVelocities;
1418
    double prevStepSize;
1419
1420
};

1421
/**
1422
 * This kernel is invoked by LangevinMiddleIntegrator to take one time step.
1423
 */
1424
class ReferenceIntegrateLangevinMiddleStepKernel : public IntegrateLangevinMiddleStepKernel {
1425
public:
1426
    ReferenceIntegrateLangevinMiddleStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateLangevinMiddleStepKernel(name, platform),
1427
1428
        data(data), dynamics(0) {
    }
1429
    ~ReferenceIntegrateLangevinMiddleStepKernel();
1430
1431
1432
1433
    /**
     * Initialize the kernel, setting up the particle masses.
     * 
     * @param system     the System this kernel will be applied to
1434
     * @param integrator the LangevinMiddleIntegrator this kernel will be used for
1435
     */
1436
    void initialize(const System& system, const LangevinMiddleIntegrator& integrator);
1437
1438
1439
1440
    /**
     * Execute the kernel.
     * 
     * @param context    the context in which to execute this kernel
1441
     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
1442
     */
1443
    void execute(ContextImpl& context, const LangevinMiddleIntegrator& integrator);
1444
1445
1446
1447
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
1448
     * @param integrator the LangevinMiddleIntegrator this kernel is being used for
1449
     */
1450
    double computeKineticEnergy(ContextImpl& context, const LangevinMiddleIntegrator& integrator);
1451
1452
private:
    ReferencePlatform::PlatformData& data;
1453
    ReferenceLangevinMiddleDynamics* dynamics;
1454
1455
1456
1457
    std::vector<double> masses;
    double prevTemp, prevFriction, prevStepSize;
};

1458
1459
1460
1461
1462
/**
 * This kernel is invoked by BrownianIntegrator to take one time step.
 */
class ReferenceIntegrateBrownianStepKernel : public IntegrateBrownianStepKernel {
public:
1463
    ReferenceIntegrateBrownianStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateBrownianStepKernel(name, platform),
1464
        data(data), dynamics(0) {
1465
    }
1466
    ~ReferenceIntegrateBrownianStepKernel();
1467
    /**
1468
     * Initialize the kernel.
1469
     * 
1470
1471
     * @param system     the System this kernel will be applied to
     * @param integrator the BrownianIntegrator this kernel will be used for
1472
     */
1473
    void initialize(const System& system, const BrownianIntegrator& integrator);
1474
1475
1476
    /**
     * Execute the kernel.
     * 
1477
1478
     * @param context    the context in which to execute this kernel
     * @param integrator the BrownianIntegrator this kernel is being used for
1479
     */
1480
    void execute(ContextImpl& context, const BrownianIntegrator& integrator);
1481
1482
1483
1484
1485
1486
1487
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the BrownianIntegrator this kernel is being used for
     */
    double computeKineticEnergy(ContextImpl& context, const BrownianIntegrator& integrator);
1488
private:
1489
    ReferencePlatform::PlatformData& data;
1490
    ReferenceBrownianDynamics* dynamics;
peastman's avatar
peastman committed
1491
    std::vector<double> masses;
1492
    double prevTemp, prevFriction, prevStepSize;
1493
1494
};

1495
1496
1497
1498
1499
1500
/**
 * 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),
1501
        data(data), dynamics(0) {
1502
1503
1504
1505
1506
1507
    }
    ~ReferenceIntegrateVariableLangevinStepKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
1508
     * @param integrator the VariableLangevinIntegrator this kernel will be used for
1509
1510
1511
1512
1513
1514
     */
    void initialize(const System& system, const VariableLangevinIntegrator& integrator);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
1515
     * @param integrator the VariableLangevinIntegrator this kernel is being used for
1516
     * @param maxTime    the maximum time beyond which the simulation should not be advanced
1517
     * @return the size of the step that was taken
1518
     */
1519
    double execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
1520
1521
1522
1523
1524
1525
1526
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the VariableLangevinIntegrator this kernel is being used for
     */
    double computeKineticEnergy(ContextImpl& context, const VariableLangevinIntegrator& integrator);
1527
1528
1529
private:
    ReferencePlatform::PlatformData& data;
    ReferenceVariableStochasticDynamics* dynamics;
peastman's avatar
peastman committed
1530
    std::vector<double> masses;
1531
1532
1533
    double prevTemp, prevFriction, prevErrorTol;
};

1534
1535
1536
1537
1538
1539
/**
 * 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),
1540
        data(data), dynamics(0) {
1541
1542
1543
1544
1545
1546
    }
    ~ReferenceIntegrateVariableVerletStepKernel();
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
1547
     * @param integrator the VariableVerletIntegrator this kernel will be used for
1548
1549
1550
1551
1552
1553
     */
    void initialize(const System& system, const VariableVerletIntegrator& integrator);
    /**
     * Execute the kernel.
     *
     * @param context    the context in which to execute this kernel
1554
     * @param integrator the VariableVerletIntegrator this kernel is being used for
1555
     * @param maxTime    the maximum time beyond which the simulation should not be advanced
1556
     * @return the size of the step that was taken
1557
     */
1558
    double execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime);
1559
1560
1561
1562
1563
1564
1565
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the VariableVerletIntegrator this kernel is being used for
     */
    double computeKineticEnergy(ContextImpl& context, const VariableVerletIntegrator& integrator);
1566
1567
1568
private:
    ReferencePlatform::PlatformData& data;
    ReferenceVariableVerletDynamics* dynamics;
peastman's avatar
peastman committed
1569
    std::vector<double> masses;
1570
    double prevErrorTol;
1571
1572
};

1573
1574
1575
1576
1577
1578
/**
 * 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),
1579
        data(data), dynamics(0) {
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
    }
    ~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);
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the CustomIntegrator this kernel is being used for
     * @param forcesAreValid if the context has been modified since the last time step, this will be
     *                       false to show that cached forces are invalid and must be recalculated.
     *                       On exit, this should specify whether the cached forces are valid at the
     *                       end of the step.
     */
    double computeKineticEnergy(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid);
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
    /**
     * 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;
peastman's avatar
peastman committed
1644
1645
    std::vector<double> masses, globalValues;
    std::vector<std::vector<OpenMM::Vec3> > perDofValues; 
1646
1647
};

Peter Eastman's avatar
Peter Eastman committed
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
/**
 * This kernel is invoked by DPDIntegrator to take one time step.
 */
class ReferenceIntegrateDPDStepKernel : public IntegrateDPDStepKernel {
public:
    ReferenceIntegrateDPDStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateDPDStepKernel(name, platform),
        data(data), dynamics(NULL) {
    }
    ~ReferenceIntegrateDPDStepKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param integrator the DPDIntegrator this kernel will be used for
     */
    void initialize(const System& system, const DPDIntegrator& integrator);
    /**
     * Execute the kernel.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the DPDIntegrator this kernel is being used for
     */
    void execute(ContextImpl& context, const DPDIntegrator& integrator);
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the DPDIntegrator this kernel is being used for
     */
    double computeKineticEnergy(ContextImpl& context, const DPDIntegrator& integrator);
private:
    ReferencePlatform::PlatformData& data;
    ReferenceDPDDynamics* dynamics;
    std::vector<double> masses;
};

1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
/**
 * This kernel is invoked by QTBIntegrator to take one time step.
 */
class ReferenceIntegrateQTBStepKernel : public IntegrateQTBStepKernel {
public:
    ReferenceIntegrateQTBStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateQTBStepKernel(name, platform),
        data(data), dynamics(NULL), hasInitialized(false) {
    }
    ~ReferenceIntegrateQTBStepKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param integrator the QTBIntegrator this kernel will be used for
     */
    void initialize(const System& system, const QTBIntegrator& integrator);
    /**
     * Execute the kernel.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the QTBIntegrator this kernel is being used for
     */
    void execute(ContextImpl& context, const QTBIntegrator& integrator);
    /**
     * Compute the kinetic energy.
     * 
     * @param context    the context in which to execute this kernel
     * @param integrator the QTBIntegrator this kernel is being used for
     */
    double computeKineticEnergy(ContextImpl& context, const QTBIntegrator& integrator);
    /**
     * Get the adapted friction coefficients for a particle.
     * 
     * @param context    the context in which to execute this kernel
     * @param particle   the index of the particle for which to get the friction
     * @param friction   the adapted friction coefficients used in generating the
     *                   random force
     */
    void getAdaptedFriction(ContextImpl& context, int particle, std::vector<double>& friction) const;
    /**
     * Set the adapted friction coefficients for a particle.  This affects the
     * specified particle, and all others that have the same type.
     * 
     * @param context    the context in which to execute this kernel
     * @param particle   the index of the particle for which to get the friction
     * @param friction   the adapted friction coefficients used in generating the
     *                   random force
     */
    void setAdaptedFriction(ContextImpl& context, int particle, const std::vector<double>& friction);
    /**
     * Write the adapted friction to a checkpoint.
     * 
     * @param context    the context in which to execute this kernel
     * @param stream     the stream to write the checkpoint to
     */
    void createCheckpoint(ContextImpl& context, std::ostream& stream) const;
    /**
     * Load the adapted friction from a checkpoint.
     * 
     * @param context    the context in which to execute this kernel
     * @param stream     the stream to read the checkpoint from
     */
    void loadCheckpoint(ContextImpl& context, std::istream& stream);
private:
    ReferencePlatform::PlatformData& data;
    ReferenceQTBDynamics* dynamics;
    std::vector<double> masses;
    bool hasInitialized;
};

1754
/**
Peter Eastman's avatar
Peter Eastman committed
1755
 * This kernel is invoked by AndersenThermostat at the start of each time step to adjust the particle velocities.
1756
1757
1758
 */
class ReferenceApplyAndersenThermostatKernel : public ApplyAndersenThermostatKernel {
public:
1759
    ReferenceApplyAndersenThermostatKernel(std::string name, const Platform& platform) : ApplyAndersenThermostatKernel(name, platform), thermostat(0) {
1760
    }
1761
    ~ReferenceApplyAndersenThermostatKernel();
1762
    /**
1763
     * Initialize the kernel.
1764
     * 
1765
1766
     * @param system     the System this kernel will be applied to
     * @param thermostat the AndersenThermostat this kernel will be used for
1767
     */
1768
    void initialize(const System& system, const AndersenThermostat& thermostat);
1769
1770
1771
    /**
     * Execute the kernel.
     * 
1772
     * @param context    the context in which to execute this kernel
1773
     */
1774
    void execute(ContextImpl& context);
1775
1776
private:
    ReferenceAndersenThermostat* thermostat;
1777
    std::vector<std::vector<int> > particleGroups;
peastman's avatar
peastman committed
1778
    std::vector<double> masses;
1779
1780
};

1781
1782
1783
1784
1785
/**
 * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
 */
class ReferenceApplyMonteCarloBarostatKernel : public ApplyMonteCarloBarostatKernel {
public:
1786
    ReferenceApplyMonteCarloBarostatKernel(std::string name, const Platform& platform) : ApplyMonteCarloBarostatKernel(name, platform), barostat(NULL) {
1787
1788
1789
1790
1791
    }
    ~ReferenceApplyMonteCarloBarostatKernel();
    /**
     * Initialize the kernel.
     *
1792
1793
     * @param system          the System this kernel will be applied to
     * @param barostat        the MonteCarloBarostat this kernel will be used for
1794
     * @param rigidMolecules  whether molecules should be kept rigid while scaling coordinates
1795
1796
     * @param components      the number of box components the barostat operates one (1 for isotropic scaling,
     *                        3 for anisotropic, 6 for both lengths and angles)
1797
     */
1798
    void initialize(const System& system, const Force& barostat, int components, bool rigidMolecules=true);
1799
1800
1801
1802
1803
1804
1805
    /**
     * Save the coordinates before attempting a Monte Carlo step.  This allows us to restore them
     * if the step is rejected.
     *
     * @param context    the context in which to execute this kernel
     */
    void saveCoordinates(ContextImpl& context);
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
    /**
     * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
     * This version scales the x, y, and z positions independently.
     * This is called BEFORE the periodic box size is modified.  It should begin by translating each particle
     * or cluster into the first periodic box, so that coordinates will still be correct after the box size
     * is changed.
     *
     * @param context    the context in which to execute this kernel
     * @param scaleX     the scale factor by which to multiply particle x-coordinate
     * @param scaleY     the scale factor by which to multiply particle y-coordinate
     * @param scaleZ     the scale factor by which to multiply particle z-coordinate
     */
1818
    void scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ);
1819
    /**
1820
1821
     * Reject the most recent Monte Carlo step, restoring the particle positions to where they were when
     * saveCoordinates() was last called.
1822
1823
1824
1825
     *
     * @param context    the context in which to execute this kernel
     */
    void restoreCoordinates(ContextImpl& context);
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
    /**
     * Compute the kinetic energy of the system.  If initialize() was called with rigidMolecules=true, this
     * should include only the translational center of mass motion of molecules.  Otherwise it should include
     * the total kinetic energy of all particles.  This is used when computing instantaneous pressure.
     * 
     * @param context    the context in which to execute this kernel
     * @param ke         a vector to store the kinetic energy components into.  On output, its length will
     *                   equal the number of components passed to initialize().
     */
    void computeKineticEnergy(ContextImpl& context, std::vector<double>& ke);
1836
private:
1837
    bool rigidMolecules;
1838
    int components;
1839
1840
1841
    ReferenceMonteCarloBarostat* barostat;
};

1842
1843
1844
1845
1846
/**
 * This kernel is invoked to remove center of mass motion from the system.
 */
class ReferenceRemoveCMMotionKernel : public RemoveCMMotionKernel {
public:
1847
    ReferenceRemoveCMMotionKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : RemoveCMMotionKernel(name, platform), data(data) {
1848
1849
    }
    /**
Peter Eastman's avatar
Peter Eastman committed
1850
     * Initialize the kernel, setting up the particle masses.
1851
     * 
1852
1853
     * @param system     the System this kernel will be applied to
     * @param force      the CMMotionRemover this kernel will be used for
1854
     */
1855
    void initialize(const System& system, const CMMotionRemover& force);
1856
1857
1858
    /**
     * Execute the kernel.
     * 
1859
     * @param context    the context in which to execute this kernel
1860
     */
1861
    void execute(ContextImpl& context);
1862
private:
1863
    ReferencePlatform::PlatformData& data;
1864
    std::vector<double> masses;
1865
    int frequency;
1866
1867
};

1868
1869
1870
1871
1872
1873
1874
1875
1876
/**
 * This kernel is invoked by ATMForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcATMForceKernel : public CalcATMForceKernel {
public:
    ReferenceCalcATMForceKernel(std::string name, const Platform& platform) : CalcATMForceKernel(name, platform) {
    }
    /**
     * Initialize the kernel.
1877
     *
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
     * @param system     the System this kernel will be applied to
     * @param force      the ATMForce this kernel will be used for
     */
    void initialize(const System& system, const ATMForce& force);
    /**
     * Scale the forces from the inner contexts and apply them to the main context.
     *
     * @param context        the context in which to execute this kernel
     * @param innerContext1  the first inner context
     * @param innerContext2  the second inner context
     * @param dEdu0          the derivative of the final energy with respect to the first inner context's energy
     * @param dEdu1          the derivative of the final energy with respect to the second inner context's energy
     * @param energyParamDerivs  derivatives of the final energy with respect to global parameters
     */
    void applyForces(ContextImpl& context, ContextImpl& innerContext0, ContextImpl& innerContext1,
                     double dEdu0, double dEdu1, const std::map<std::string, double>& energyParamDerivs);
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the ATMForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const ATMForce& force);
    /**
     * Copy state information to the inner contexts.
     *
     * @param context        the context in which to execute this kernel
     * @param innerContext1  the first context created by the ATMForce for computing displaced energy
     * @param innerContext2  the second context created by the ATMForce for computing displaced energy
     */
    void copyState(ContextImpl& context, ContextImpl& innerContext0, ContextImpl& innerContext1);
private:
    int numParticles;
1911
1912
1913
1914
1915
1916
    std::vector<Vec3> displ1, displ0;
    std::vector<Vec3> displacement1, displacement0;
    std::vector<int> pj1, pi1, pj0, pi0;
    void setDisplacements(std::vector<Vec3>& pos);
    void displForces(std::vector<Vec3>& force0, std::vector<Vec3>& force1);
    void loadParams(int numParticles, const ATMForce& force);
1917
1918
};

1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
/**
 * This kernel is invoked by CustomCPPForceImpl to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcCustomCPPForceKernel : public CalcCustomCPPForceKernel {
public:
    ReferenceCalcCustomCPPForceKernel(std::string name, const Platform& platform) : CalcCustomCPPForceKernel(name, platform), force(NULL) {
    }
    /**
     * Initialize the kernel.
     *
     * @param system     the System this kernel will be applied to
     * @param force      the CustomCPPForceImpl this kernel will be used for
     */
    void initialize(const System& system, CustomCPPForceImpl& 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:
    CustomCPPForceImpl* force;
    std::vector<Vec3> forces;
};

1947
1948
1949
} // namespace OpenMM

#endif /*OPENMM_REFERENCEKERNELS_H_*/