AmoebaReferenceKernels.h 29.8 KB
Newer Older
1
2
3
4
#ifndef AMOEBA_OPENMM_REFERENCE_KERNELS_H_
#define AMOEBA_OPENMM_REFERENCE_KERNELS_H_

/* -------------------------------------------------------------------------- *
5
 *                              OpenMMAmoeba                                  *
6
7
8
9
10
11
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
peastman's avatar
peastman committed
12
 * Portions copyright (c) 2008-2018 Stanford University and the Authors.      *
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
 * Authors:                                                                   *
 * Contributors:                                                              *
 *                                                                            *
 * This program is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU Lesser General Public License as published   *
 * by the Free Software Foundation, either version 3 of the License, or       *
 * (at your option) any later version.                                        *
 *                                                                            *
 * This program is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
 * GNU Lesser General Public License for more details.                        *
 *                                                                            *
 * You should have received a copy of the GNU Lesser General Public License   *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
 * -------------------------------------------------------------------------- */

#include "openmm/System.h"
31
#include "openmm/amoebaKernels.h"
32
#include "openmm/AmoebaMultipoleForce.h"
peastman's avatar
peastman committed
33
#include "openmm/HippoNonbondedForce.h"
34
#include "AmoebaReferenceMultipoleForce.h"
peastman's avatar
peastman committed
35
#include "AmoebaReferenceHippoNonbondedForce.h"
36
37
#include "ReferenceNeighborList.h"
#include "SimTKOpenMMRealType.h"
38
39
40
41

namespace OpenMM {

/**
42
 * This kernel is invoked by AmoebaBondForce to calculate the forces acting on the system and the energy of the system.
43
 */
44
class ReferenceCalcAmoebaBondForceKernel : public CalcAmoebaBondForceKernel {
45
public:
46
    ReferenceCalcAmoebaBondForceKernel(std::string name, 
47
                                               const Platform& platform,
48
                                               const System& system);
49
    ~ReferenceCalcAmoebaBondForceKernel();
50
51
52
53
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
54
     * @param force      the AmoebaBondForce this kernel will be used for
55
     */
56
    void initialize(const System& system, const AmoebaBondForce& force);
57
58
59
60
61
62
63
64
65
    /**
     * 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);
66
67
68
69
70
71
72
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the AmoebaBondForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const AmoebaBondForce& force);
73
74
75
76
private:
    int numBonds;
    std::vector<int>   particle1;
    std::vector<int>   particle2;
peastman's avatar
peastman committed
77
78
79
80
    std::vector<double> length;
    std::vector<double> kQuadratic;
    double globalBondCubic;
    double globalBondQuartic;
81
    const System& system;
82
    bool usePeriodic;
83
84
85
};

/**
86
 * This kernel is invoked by AmoebaAngleForce to calculate the forces acting on the system and the energy of the system.
87
 */
88
class ReferenceCalcAmoebaAngleForceKernel : public CalcAmoebaAngleForceKernel {
89
public:
90
    ReferenceCalcAmoebaAngleForceKernel(std::string name, const Platform& platform, const System& system);
91
    ~ReferenceCalcAmoebaAngleForceKernel();
92
93
94
95
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
96
     * @param force      the AmoebaAngleForce this kernel will be used for
97
     */
98
    void initialize(const System& system, const AmoebaAngleForce& force);
99
100
101
102
103
104
105
106
107
    /**
     * 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);
108
109
110
111
112
113
114
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the AmoebaAngleForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const AmoebaAngleForce& force);
115
116
117
118
119
private:
    int numAngles;
    std::vector<int>   particle1;
    std::vector<int>   particle2;
    std::vector<int>   particle3;
peastman's avatar
peastman committed
120
121
122
123
124
125
    std::vector<double> angle;
    std::vector<double> kQuadratic;
    double globalAngleCubic;
    double globalAngleQuartic;
    double globalAnglePentic;
    double globalAngleSextic;
126
    const System& system;
127
    bool usePeriodic;
128
129
};

130
/**
131
 * This kernel is invoked by AmoebaInPlaneAngleForce to calculate the forces acting on the system and the energy of the system.
132
 */
133
class ReferenceCalcAmoebaInPlaneAngleForceKernel : public CalcAmoebaInPlaneAngleForceKernel {
134
public:
135
    ReferenceCalcAmoebaInPlaneAngleForceKernel(std::string name, const Platform& platform, const System& system);
136
    ~ReferenceCalcAmoebaInPlaneAngleForceKernel();
137
138
139
140
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
141
     * @param force      the AmoebaInPlaneAngleForce this kernel will be used for
142
     */
143
    void initialize(const System& system, const AmoebaInPlaneAngleForce& force);
144
145
146
147
148
149
150
151
152
    /**
     * 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);
153
154
155
156
157
158
159
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the AmoebaInPlaneAngleForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const AmoebaInPlaneAngleForce& force);
160
161
162
163
164
165
private:
    int numAngles;
    std::vector<int>   particle1;
    std::vector<int>   particle2;
    std::vector<int>   particle3;
    std::vector<int>   particle4;
peastman's avatar
peastman committed
166
167
168
169
170
171
    std::vector<double> angle;
    std::vector<double> kQuadratic;
    double globalInPlaneAngleCubic;
    double globalInPlaneAngleQuartic;
    double globalInPlaneAnglePentic;
    double globalInPlaneAngleSextic;
172
    const System& system;
173
    bool usePeriodic;
174
175
};

176
177
178
179
180
/**
 * This kernel is invoked by AmoebaPiTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcAmoebaPiTorsionForceKernel : public CalcAmoebaPiTorsionForceKernel {
public:
181
    ReferenceCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, const System& system);
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
    ~ReferenceCalcAmoebaPiTorsionForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param force      the AmoebaPiTorsionForce this kernel will be used for
     */
    void initialize(const System& system, const AmoebaPiTorsionForce& 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);
199
200
201
202
203
204
205
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the AmoebaPiTorsionForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const AmoebaPiTorsionForce& force);
206
207
208
209
210
211
212
213
private:
    int numPiTorsions;
    std::vector<int>   particle1;
    std::vector<int>   particle2;
    std::vector<int>   particle3;
    std::vector<int>   particle4;
    std::vector<int>   particle5;
    std::vector<int>   particle6;
peastman's avatar
peastman committed
214
    std::vector<double> kTorsion;
215
    const System& system;
216
    bool usePeriodic;
217
218
};

219
220
221
222
223
/**
 * This kernel is invoked by AmoebaStretchBendForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcAmoebaStretchBendForceKernel : public CalcAmoebaStretchBendForceKernel {
public:
224
    ReferenceCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, const System& system);
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
    ~ReferenceCalcAmoebaStretchBendForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param force      the AmoebaStretchBendForce this kernel will be used for
     */
    void initialize(const System& system, const AmoebaStretchBendForce& 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);
242
243
244
245
246
247
248
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the AmoebaStretchBendForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const AmoebaStretchBendForce& force);
249
250
251
252
253
private:
    int numStretchBends;
    std::vector<int>   particle1;
    std::vector<int>   particle2;
    std::vector<int>   particle3;
peastman's avatar
peastman committed
254
255
256
257
258
    std::vector<double> lengthABParameters;
    std::vector<double> lengthCBParameters;
    std::vector<double> angleParameters;
    std::vector<double> k1Parameters;
    std::vector<double> k2Parameters;
259
    const System& system;
260
    bool usePeriodic;
261
};
262
263
264
265
266
267

/**
 * This kernel is invoked by AmoebaOutOfPlaneBendForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcAmoebaOutOfPlaneBendForceKernel : public CalcAmoebaOutOfPlaneBendForceKernel {
public:
268
    ReferenceCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, const System& system);
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
    ~ReferenceCalcAmoebaOutOfPlaneBendForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param force      the AmoebaOutOfPlaneBendForce this kernel will be used for
     */
    void initialize(const System& system, const AmoebaOutOfPlaneBendForce& 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);
286
287
288
289
290
291
292
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the AmoebaOutOfPlaneBendForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const AmoebaOutOfPlaneBendForce& force);
293
294
295
296
297
298
private:
    int numOutOfPlaneBends;
    std::vector<int>   particle1;
    std::vector<int>   particle2;
    std::vector<int>   particle3;
    std::vector<int>   particle4;
peastman's avatar
peastman committed
299
300
301
302
303
    std::vector<double> kParameters;
    double globalOutOfPlaneBendAngleCubic;
    double globalOutOfPlaneBendAngleQuartic;
    double globalOutOfPlaneBendAnglePentic;
    double globalOutOfPlaneBendAngleSextic;
304
    const System& system;
305
    bool usePeriodic;
306
307
};

308
309
310
311
312
/**
 * This kernel is invoked by AmoebaTorsionTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcAmoebaTorsionTorsionForceKernel : public CalcAmoebaTorsionTorsionForceKernel {
public:
313
    ReferenceCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, const System& system);
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
    ~ReferenceCalcAmoebaTorsionTorsionForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param force      the AmoebaTorsionTorsionForce this kernel will be used for
     */
    void initialize(const System& system, const AmoebaTorsionTorsionForce& force);
    /**
     * Execute the kernel to calculate the forces and/or energy.
     *
     * @param context        the context in which to execute this kernel
     * @param includeForces  true if forces should be calculated
     * @param includeEnergy  true if the energy should be calculated
     * @return the potential energy due to the force
     */
    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
private:
    int numTorsionTorsions;
    std::vector<int>   particle1;
    std::vector<int>   particle2;
    std::vector<int>   particle3;
    std::vector<int>   particle4;
    std::vector<int>   particle5;
    std::vector<int>   chiralCheckAtom;
    std::vector<int>   gridIndices;

    int numTorsionTorsionGrids;
peastman's avatar
peastman committed
342
    std::vector< std::vector< std::vector< std::vector<double> > > > torsionTorsionGrids;
343

344
    const System& system;
345
    bool usePeriodic;
346
347
};

348
349
350
351
352
/**
 * This kernel is invoked by AmoebaMultipoleForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcAmoebaMultipoleForceKernel : public CalcAmoebaMultipoleForceKernel {
public:
353
    ReferenceCalcAmoebaMultipoleForceKernel(std::string name, const Platform& platform, const System& system);
354
355
356
357
358
359
360
361
    ~ReferenceCalcAmoebaMultipoleForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param force      the AmoebaMultipoleForce this kernel will be used for
     */
    void initialize(const System& system, const AmoebaMultipoleForce& force);
362
363
364
365
366
367
368
    /**
     * Setup for AmoebaReferenceMultipoleForce instance. 
     *
     * @param context        the current context
     *
     * @return pointer to initialized instance of AmoebaReferenceMultipoleForce
     */
369
    AmoebaReferenceMultipoleForce* setupAmoebaReferenceMultipoleForce(ContextImpl& context);
370
371
372
373
374
375
376
377
378
    /**
     * 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);
379
380
381
382
383
384
385
    /**
     * Get the induced dipole moments of all particles.
     * 
     * @param context    the Context for which to get the induced dipoles
     * @param dipoles    the induced dipole moment of particle i is stored into the i'th element
     */
    void getInducedDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
386
    /**
387
388
389
390
391
     * Get the fixed dipole moments of all particles in the global reference frame.
     * 
     * @param context    the Context for which to get the fixed dipoles
     * @param dipoles    the fixed dipole moment of particle i is stored into the i'th element
     */
392
    void getLabFramePermanentDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
393
394
395
396
397
398
399
    /**
     * Get the total dipole moments of all particles in the global reference frame.
     * 
     * @param context    the Context for which to get the fixed dipoles
     * @param dipoles    the fixed dipole moment of particle i is stored into the i'th element
     */
    void getTotalDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
400
    /** 
401
     * Calculate the electrostatic potential given vector of grid coordinates.
402
     *
403
404
     * @param context                      context
     * @param inputGrid                    input grid coordinates
405
406
407
     * @param outputElectrostaticPotential output potential 
     */
    void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
408
                                   std::vector< double >& outputElectrostaticPotential);
409
410

    /**
411
     * Get the system multipole moments.
412
     *
413
     * @param context                context 
414
     * @param outputMultipoleMoments vector of multipole moments:
415
                                     (charge,
416
417
418
                                      dipole_x, dipole_y, dipole_z,
                                      quadrupole_xx, quadrupole_xy, quadrupole_xz,
                                      quadrupole_yx, quadrupole_yy, quadrupole_yz,
419
                                      quadrupole_zx, quadrupole_zy, quadrupole_zz)
420
     */
Lee-Ping Wang's avatar
Lee-Ping Wang committed
421
    void getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMoments);
422
423
424
425
426
427
428
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the AmoebaMultipoleForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const AmoebaMultipoleForce& force);
429
430
431
432
433
434
435
436
437
    /**
     * 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;
438

439
private:
440

441
    int numMultipoles;
442
    AmoebaMultipoleForce::NonbondedMethod nonbondedMethod;
443
    AmoebaMultipoleForce::PolarizationType polarizationType;
peastman's avatar
peastman committed
444
445
446
447
448
449
    std::vector<double> charges;
    std::vector<double> dipoles;
    std::vector<double> quadrupoles;
    std::vector<double> tholes;
    std::vector<double> dampingFactors;
    std::vector<double> polarity;
450
    std::vector<int>   axisTypes;
451
452
453
    std::vector<int>   multipoleAtomZs;
    std::vector<int>   multipoleAtomXs;
    std::vector<int>   multipoleAtomYs;
454
455
456
    std::vector< std::vector< std::vector<int> > > multipoleAtomCovalentInfo;

    int mutualInducedMaxIterations;
peastman's avatar
peastman committed
457
    double mutualInducedTargetEpsilon;
458
    std::vector<double> extrapolationCoefficients;
459

460
    bool usePme;
peastman's avatar
peastman committed
461
462
    double alphaEwald;
    double cutoffDistance;
463
464
    std::vector<int> pmeGridDimension;

465
    const System& system;
466
467
};

Mark Friedrichs's avatar
Mark Friedrichs committed
468
469
470
471
472
/**
 * This kernel is invoked to calculate the vdw forces acting on the system and the energy of the system.
 */
class ReferenceCalcAmoebaVdwForceKernel : public CalcAmoebaVdwForceKernel {
public:
473
    ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, const System& system);
Mark Friedrichs's avatar
Mark Friedrichs committed
474
475
476
477
478
    ~ReferenceCalcAmoebaVdwForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
479
     * @param force      the AmoebaVdwForce this kernel will be used for
Mark Friedrichs's avatar
Mark Friedrichs committed
480
481
482
483
484
485
486
487
488
489
490
     */
    void initialize(const System& system, const AmoebaVdwForce& 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);
491
492
493
494
495
496
497
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the AmoebaVdwForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const AmoebaVdwForce& force);
Mark Friedrichs's avatar
Mark Friedrichs committed
498
499
private:
    int numParticles;
500
    int useCutoff;
501
502
    int usePBC;
    double cutoff;
503
    double dispersionCoefficient;
Mark Friedrichs's avatar
Mark Friedrichs committed
504
    std::vector<int> indexIVs;
505
    std::vector< std::set<int> > allExclusions;
peastman's avatar
peastman committed
506
507
508
    std::vector<double> sigmas;
    std::vector<double> epsilons;
    std::vector<double> reductions;
Mark Friedrichs's avatar
Mark Friedrichs committed
509
510
    std::string sigmaCombiningRule;
    std::string epsilonCombiningRule;
511
    const System& system;
512
    NeighborList* neighborList;
Mark Friedrichs's avatar
Mark Friedrichs committed
513
514
};

Mark Friedrichs's avatar
Mark Friedrichs committed
515
516
517
518
519
/**
 * This kernel is invoked to calculate the WCA dispersion forces acting on the system and the energy of the system.
 */
class ReferenceCalcAmoebaWcaDispersionForceKernel : public CalcAmoebaWcaDispersionForceKernel {
public:
520
    ReferenceCalcAmoebaWcaDispersionForceKernel(std::string name, const Platform& platform, const System& system);
Mark Friedrichs's avatar
Mark Friedrichs committed
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
    ~ReferenceCalcAmoebaWcaDispersionForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param force      the AmoebaMultipoleForce this kernel will be used for
     */
    void initialize(const System& system, const AmoebaWcaDispersionForce& 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);
538
539
540
541
542
543
544
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the AmoebaWcaDispersionForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const AmoebaWcaDispersionForce& force);
Mark Friedrichs's avatar
Mark Friedrichs committed
545
546
547
private:

    int numParticles;
peastman's avatar
peastman committed
548
549
550
551
552
553
554
555
556
557
558
    std::vector<double> radii;
    std::vector<double> epsilons;
    double epso; 
    double epsh; 
    double rmino; 
    double rminh; 
    double awater; 
    double shctd; 
    double dispoff;
    double slevy;
    double totalMaximumDispersionEnergy;
559
    const System& system;
Mark Friedrichs's avatar
Mark Friedrichs committed
560
};
561

562
563
564
565
566
/**
 * This kernel is invoked to calculate the Gerneralized Kirkwood forces acting on the system and the energy of the system.
 */
class ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel : public CalcAmoebaGeneralizedKirkwoodForceKernel {
public:
567
    ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, const System& system);
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
    ~ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param force      the AmoebaMultipoleForce this kernel will be used for
     */
    void initialize(const System& system, const AmoebaGeneralizedKirkwoodForce& 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);
 
    /**
587
     *  Get the 'include cavity term' flag.
588
589
590
     *
     *  @return includeCavityTerm
     */
591
    int getIncludeCavityTerm() const;
592
593

    /**
594
     *  Get the number of particles.
595
596
597
     *
     *  @return number of particles
     */
598
    int getNumParticles() const;
599
600

    /**
601
     *  Get Direct Polarization flag.
602
603
604
605
     *
     *  @return directPolarization
     *
     */
606
    int getDirectPolarization() const;
607
608

    /**
609
     *  Get the solute dielectric.
610
611
612
613
     *
     *  @return soluteDielectric
     *
     */
peastman's avatar
peastman committed
614
    double getSoluteDielectric() const;
615
616

    /**
617
     *  Get the solvent dielectric.
618
619
620
621
     *
     *  @return solventDielectric
     *
     */
peastman's avatar
peastman committed
622
    double getSolventDielectric() const;
623
624

    /**
625
     *  Get the dielectric offset.
626
627
628
629
     *
     *  @return dielectricOffset
     *
     */
peastman's avatar
peastman committed
630
    double getDielectricOffset() const;
631
632

    /**
633
     *  Get the probe radius.
634
635
636
637
     *
     *  @return probeRadius
     *
     */
peastman's avatar
peastman committed
638
    double getProbeRadius() const;
639
640

    /**
641
     *  Get the surface area factor.
642
643
644
645
     *
     *  @return surfaceAreaFactor
     *
     */
peastman's avatar
peastman committed
646
    double getSurfaceAreaFactor() const;
647
648

    /**
649
     *  Get the vector of particle radii.
650
651
652
653
     *
     *  @param atomicRadii vector of atomic radii
     *
     */
peastman's avatar
peastman committed
654
    void getAtomicRadii(std::vector<double>& atomicRadii) const;
655
656

    /**
657
     *  Get the vector of scale factors.
658
659
660
661
     *
     *  @param scaleFactors vector of scale factors
     *
     */
peastman's avatar
peastman committed
662
    void getScaleFactors(std::vector<double>& scaleFactors) const;
663
664

    /**
665
     *  Get the vector of charges.
666
667
668
669
     *
     *  @param charges vector of charges
     *
     */
peastman's avatar
peastman committed
670
    void getCharges(std::vector<double>& charges) const;
671

672
673
674
675
676
677
678
679
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the AmoebaGeneralizedKirkwoodForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const AmoebaGeneralizedKirkwoodForce& force);

680
681
682
private:

    int numParticles;
peastman's avatar
peastman committed
683
684
685
686
687
688
689
690
    std::vector<double> atomicRadii;
    std::vector<double> scaleFactors;
    std::vector<double> charges;
    double soluteDielectric;
    double solventDielectric;
    double dielectricOffset;
    double probeRadius;
    double surfaceAreaFactor;
691
692
    int includeCavityTerm;
    int directPolarization;
693
    const System& system;
694
695
};

peastman's avatar
peastman committed
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
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
/**
 * This kernel is invoked by HippoNonbondedForce to calculate the forces acting on the system and the energy of the system.
 */
class ReferenceCalcHippoNonbondedForceKernel : public CalcHippoNonbondedForceKernel {
public:
    ReferenceCalcHippoNonbondedForceKernel(std::string name, const Platform& platform, const System& system);
    ~ReferenceCalcHippoNonbondedForceKernel();
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
     * @param force      the HippoNonbondedForce this kernel will be used for
     */
    void initialize(const System& system, const HippoNonbondedForce& force);
    /**
     * Setup for AmoebaReferenceHippoNonbondedForce instance. 
     *
     * @param context        the current context
     *
     * @return pointer to initialized instance of AmoebaReferenceHippoNonbondedForce
     */
    void setupAmoebaReferenceHippoNonbondedForce(ContextImpl& context);
    /**
     * 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);
    /**
     * Get the induced dipole moments of all particles.
     * 
     * @param context    the Context for which to get the induced dipoles
     * @param dipoles    the induced dipole moment of particle i is stored into the i'th element
     */
    void getInducedDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
    /**
     * Get the fixed dipole moments of all particles in the global reference frame.
     * 
     * @param context    the Context for which to get the fixed dipoles
     * @param dipoles    the fixed dipole moment of particle i is stored into the i'th element
     */
    void getLabFramePermanentDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
    /**
     * Get the total dipole moments of all particles in the global reference frame.
     * 
     * @param context    the Context for which to get the fixed dipoles
     * @param dipoles    the fixed dipole moment of particle i is stored into the i'th element
     */
    void getTotalDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
    /**
     * Copy changed parameters over to a context.
     *
     * @param context    the context to copy parameters to
     * @param force      the HippoNonbondedForce to copy the parameters from
     */
    void copyParametersToContext(ContextImpl& context, const HippoNonbondedForce& force);
    /**
     * 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 parameters being used for dispersion 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 getDPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;

private:

    AmoebaReferenceHippoNonbondedForce* ixn;
    int numParticles;
};

780
781
782
} // namespace OpenMM

#endif /*AMOEBA_OPENMM_REFERENCE_KERNELS_H*/