AmoebaCudaKernels.h 26.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#ifndef AMOEBA_OPENMM_CUDAKERNELS_H_
#define AMOEBA_OPENMM_CUDAKERNELS_H_

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

#include "openmm/amoebaKernels.h"
#include "openmm/kernels.h"
#include "openmm/System.h"
#include "CudaArray.h"
34
35
36
#include "CudaContext.h"
#include "CudaSort.h"
#include <cufft.h>
37
38
39

namespace OpenMM {

40
41
class CudaCalcAmoebaGeneralizedKirkwoodForceKernel;

42
/**
43
 * This kernel is invoked by AmoebaBondForce to calculate the forces acting on the system and the energy of the system.
44
 */
45
class CudaCalcAmoebaBondForceKernel : public CalcAmoebaBondForceKernel {
46
public:
47
    CudaCalcAmoebaBondForceKernel(std::string name, 
48
49
                                          const Platform& platform,
                                          CudaContext& cu,
50
                                          const System& system);
51
52
53
54
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
55
     * @param force      the AmoebaBondForce this kernel will be used for
56
     */
57
    void initialize(const System& system, const AmoebaBondForce& force);
58
59
60
61
62
63
64
65
66
    /**
     * 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);
67
68
69
70
71
72
73
    /**
     * 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);
74
75
76
77
private:
    class ForceInfo;
    int numBonds;
    CudaContext& cu;
78
    const System& system;
79
    CudaArray params;
80
81
82
};

/**
83
 * This kernel is invoked by AmoebaAngleForce to calculate the forces acting on the system and the energy of the system.
84
 */
85
class CudaCalcAmoebaAngleForceKernel : public CalcAmoebaAngleForceKernel {
86
public:
87
    CudaCalcAmoebaAngleForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system);
88
89
90
91
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
92
     * @param force      the AmoebaAngleForce this kernel will be used for
93
     */
94
    void initialize(const System& system, const AmoebaAngleForce& force);
95
96
97
98
99
100
101
102
103
    /**
     * 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);
104
105
106
107
108
109
110
    /**
     * 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);
111
112
113
114
private:
    class ForceInfo;
    int numAngles;
    CudaContext& cu;
115
    const System& system;
116
    CudaArray params;
117
118
119
};

/**
120
 * This kernel is invoked by AmoebaInPlaneAngleForce to calculate the forces acting on the system and the energy of the system.
121
 */
122
class CudaCalcAmoebaInPlaneAngleForceKernel : public CalcAmoebaInPlaneAngleForceKernel {
123
public:
124
    CudaCalcAmoebaInPlaneAngleForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system);
125
126
127
128
    /**
     * Initialize the kernel.
     * 
     * @param system     the System this kernel will be applied to
129
     * @param force      the AmoebaInPlaneAngleForce this kernel will be used for
130
     */
131
    void initialize(const System& system, const AmoebaInPlaneAngleForce& force);
132
133
134
135
136
137
138
139
140
    /**
     * 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);
141
142
143
144
145
146
147
    /**
     * 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);
148
149
150
151
private:
    class ForceInfo;
    int numAngles;
    CudaContext& cu;
152
    const System& system;
153
    CudaArray params;
154
155
156
157
158
159
160
};

/**
 * This kernel is invoked by AmoebaPiTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class CudaCalcAmoebaPiTorsionForceKernel : public CalcAmoebaPiTorsionForceKernel {
public:
161
    CudaCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system);
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
    /**
     * 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);
178
179
180
181
182
183
184
    /**
     * 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);
185
186
187
188
private:
    class ForceInfo;
    int numPiTorsions;
    CudaContext& cu;
189
    const System& system;
190
    CudaArray params;
191
192
193
194
195
196
197
};

/**
 * This kernel is invoked by AmoebaStretchBendForce to calculate the forces acting on the system and the energy of the system.
 */
class CudaCalcAmoebaStretchBendForceKernel : public CalcAmoebaStretchBendForceKernel {
public:
198
    CudaCalcAmoebaStretchBendForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system);
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
    /**
     * 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);
215
216
217
218
219
220
221
    /**
     * 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);
222
223
224
225
private:
    class ForceInfo;
    int numStretchBends;
    CudaContext& cu;
226
    const System& system;
227
228
    CudaArray params1; // Equilibrium values
    CudaArray params2; // force constants
229
230
231
232
233
234
235
};

/**
 * This kernel is invoked by AmoebaOutOfPlaneBendForce to calculate the forces acting on the system and the energy of the system.
 */
class CudaCalcAmoebaOutOfPlaneBendForceKernel : public CalcAmoebaOutOfPlaneBendForceKernel {
public:
236
    CudaCalcAmoebaOutOfPlaneBendForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system);
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
    /**
     * 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);
253
254
255
256
257
258
259
    /**
     * 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);
260
261
262
263
private:
    class ForceInfo;
    int numOutOfPlaneBends;
    CudaContext& cu;
264
    const System& system;
265
    CudaArray params;
266
267
268
269
270
271
272
};

/**
 * This kernel is invoked by AmoebaTorsionTorsionForce to calculate the forces acting on the system and the energy of the system.
 */
class CudaCalcAmoebaTorsionTorsionForceKernel : public CalcAmoebaTorsionTorsionForceKernel {
public:
273
    CudaCalcAmoebaTorsionTorsionForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system);
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
    /**
     * 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:
    class ForceInfo;
    int numTorsionTorsions;
    int numTorsionTorsionGrids;
    CudaContext& cu;
295
    const System& system;
296
297
298
    CudaArray gridValues;
    CudaArray gridParams;
    CudaArray torsionParams;
299
300
301
302
303
304
305
};

/**
 * This kernel is invoked by AmoebaMultipoleForce to calculate the forces acting on the system and the energy of the system.
 */
class CudaCalcAmoebaMultipoleForceKernel : public CalcAmoebaMultipoleForceKernel {
public:
306
    CudaCalcAmoebaMultipoleForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system);
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
    ~CudaCalcAmoebaMultipoleForceKernel();
    /**
     * 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);
    /**
     * 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);
324
325
326
327
328
329
     /**
     * Get the LabFrame 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
     */
330
    void getLabFramePermanentDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
331
332
333
334
335
336
337
    /**
     * 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);
338
339
340
341
342
343
344
    /**
     * Get the total 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 getTotalDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
345
346
347
348
349
350
351
352
    /**
     * Execute the kernel to calculate the electrostatic potential
     *
     * @param context        the context in which to execute this kernel
     * @param inputGrid      input grid coordinates
     * @param outputElectrostaticPotential output potential 
     */
    void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
353
                                   std::vector< double >& outputElectrostaticPotential);
354
355
356
357
358

   /** 
     * Get the system multipole moments
     *
     * @param context      context
359
360
361
362
     * @param outputMultipoleMoments (charge,
     *                                dipole_x, dipole_y, dipole_z,
     *                                quadrupole_xx, quadrupole_xy, quadrupole_xz,
     *                                quadrupole_yx, quadrupole_yy, quadrupole_yz,
363
     *                                quadrupole_zx, quadrupole_zy, quadrupole_zz)
364
     */
Lee-Ping Wang's avatar
Lee-Ping Wang committed
365
    void getSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments);
366
367
368
369
370
371
372
    /**
     * 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);
373
374
375
376
377
378
379
380
381
    /**
     * 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;
382
383
private:
    class ForceInfo;
384
385
386
387
388
    class SortTrait : public CudaSort::SortTrait {
        int getDataSize() const {return 8;}
        int getKeySize() const {return 4;}
        const char* getDataType() const {return "int2";}
        const char* getKeyType() const {return "int";}
Robert McGibbon's avatar
Robert McGibbon committed
389
390
391
        const char* getMinKey() const {return "(-2147483647 - 1)";}
        const char* getMaxKey() const {return "2147483647";}
        const char* getMaxValue() const {return "make_int2(2147483647, 2147483647)";}
392
393
        const char* getSortKey() const {return "value.y";}
    };
394
    void initializeScaleFactors();
395
    void computeInducedField(void** recipBoxVectorPointer);
peastman's avatar
peastman committed
396
    bool iterateDipolesByDIIS(int iteration);
397
    void computeExtrapolatedDipoles(void** recipBoxVectorPointer);
398
    void ensureMultipolesValid(ContextImpl& context);
Lee-Ping Wang's avatar
Lee-Ping Wang committed
399
    template <class T, class T4, class M4> void computeSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments);
400
    int numMultipoles, maxInducedIterations, maxExtrapolationOrder;
401
    int fixedFieldThreads, inducedFieldThreads, electrostaticsThreads;
402
403
    int gridSizeX, gridSizeY, gridSizeZ;
    double alpha, inducedEpsilon;
Peter Eastman's avatar
Peter Eastman committed
404
    bool usePME, hasQuadrupoles, hasInitializedScaleFactors, hasInitializedFFT, multipolesAreValid, hasCreatedEvent;
405
    AmoebaMultipoleForce::PolarizationType polarizationType;
406
    CudaContext& cu;
407
    const System& system;
408
409
    std::vector<int3> covalentFlagValues;
    std::vector<int2> polarizationFlagValues;
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
    CudaArray multipoleParticles;
    CudaArray molecularDipoles;
    CudaArray molecularQuadrupoles;
    CudaArray labFrameDipoles;
    CudaArray labFrameQuadrupoles;
    CudaArray sphericalDipoles;
    CudaArray sphericalQuadrupoles;
    CudaArray fracDipoles;
    CudaArray fracQuadrupoles;
    CudaArray field;
    CudaArray fieldPolar;
    CudaArray inducedField;
    CudaArray inducedFieldPolar;
    CudaArray torque;
    CudaArray dampingAndThole;
    CudaArray inducedDipole;
    CudaArray inducedDipolePolar;
    CudaArray inducedDipoleErrors;
    CudaArray prevDipoles;
    CudaArray prevDipolesPolar;
    CudaArray prevDipolesGk;
    CudaArray prevDipolesGkPolar;
    CudaArray prevErrors;
    CudaArray diisMatrix;
    CudaArray diisCoefficients;
    CudaArray extrapolatedDipole;
    CudaArray extrapolatedDipolePolar;
    CudaArray extrapolatedDipoleGk;
    CudaArray extrapolatedDipoleGkPolar;
    CudaArray inducedDipoleFieldGradient;
    CudaArray inducedDipoleFieldGradientPolar;
    CudaArray inducedDipoleFieldGradientGk;
    CudaArray inducedDipoleFieldGradientGkPolar;
    CudaArray extrapolatedDipoleFieldGradient;
    CudaArray extrapolatedDipoleFieldGradientPolar;
    CudaArray extrapolatedDipoleFieldGradientGk;
    CudaArray extrapolatedDipoleFieldGradientGkPolar;
    CudaArray polarizability;
    CudaArray covalentFlags;
    CudaArray polarizationGroupFlags;
    CudaArray pmeGrid;
    CudaArray pmeBsplineModuliX;
    CudaArray pmeBsplineModuliY;
    CudaArray pmeBsplineModuliZ;
    CudaArray pmeIgrid;
    CudaArray pmePhi;
    CudaArray pmePhid;
    CudaArray pmePhip;
    CudaArray pmePhidp;
    CudaArray pmeCphi;
    CudaArray lastPositions;
461
462
    CudaSort* sort;
    cufftHandle fft;
463
    CUfunction computeMomentsKernel, recordInducedDipolesKernel, computeFixedFieldKernel, computeInducedFieldKernel, updateInducedFieldKernel, electrostaticsKernel, mapTorqueKernel;
Peter Eastman's avatar
Peter Eastman committed
464
    CUfunction pmeSpreadFixedMultipolesKernel, pmeSpreadInducedDipolesKernel, pmeFinishSpreadChargeKernel, pmeConvolutionKernel;
465
    CUfunction pmeFixedPotentialKernel, pmeInducedPotentialKernel, pmeFixedForceKernel, pmeInducedForceKernel, pmeRecordInducedFieldDipolesKernel, computePotentialKernel;
Peter Eastman's avatar
Peter Eastman committed
466
    CUfunction recordDIISDipolesKernel, buildMatrixKernel, solveMatrixKernel;
467
    CUfunction initExtrapolatedKernel, iterateExtrapolatedKernel, computeExtrapolatedKernel, addExtrapolatedGradientKernel;
468
    CUfunction pmeTransformMultipolesKernel, pmeTransformPotentialKernel;
Peter Eastman's avatar
Peter Eastman committed
469
    CUevent syncEvent;
470
    CudaCalcAmoebaGeneralizedKirkwoodForceKernel* gkKernel;
471
    static const int PmeOrder = 5;
peastman's avatar
peastman committed
472
    static const int MaxPrevDIISDipoles = 20;
473
474
475
476
477
478
479
};

/**
 * This kernel is invoked by AmoebaMultipoleForce to calculate the forces acting on the system and the energy of the system.
 */
class CudaCalcAmoebaGeneralizedKirkwoodForceKernel : public CalcAmoebaGeneralizedKirkwoodForceKernel {
public:
480
    CudaCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system);
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
    /**
     * 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);
497
498
499
500
501
502
503
504
    /**
     * Perform the computation of Born radii.
     */
    void computeBornRadii();
    /**
     * Perform the final parts of the force/energy computation.
     */
    void finishComputation(CudaArray& torque, CudaArray& labFrameDipoles, CudaArray& labFrameQuadrupoles, CudaArray& inducedDipole, CudaArray& inducedDipolePolar, CudaArray& dampingAndThole, CudaArray& covalentFlags, CudaArray& polarizationGroupFlags);
505
    CudaArray& getBornRadii() {
506
507
        return bornRadii;
    }
508
    CudaArray& getField() {
509
510
        return field;
    }
511
    CudaArray& getInducedField() {
512
513
        return inducedField;
    }
514
    CudaArray& getInducedFieldPolar() {
515
516
        return inducedFieldPolar;
    }
517
    CudaArray& getInducedDipoles() {
518
519
        return inducedDipoleS;
    }
520
    CudaArray& getInducedDipolesPolar() {
521
522
        return inducedDipolePolarS;
    }
523
524
525
526
527
528
529
    /**
     * 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);
530
531
532
private:
    class ForceInfo;
    CudaContext& cu;
533
    const System& system;
534
    bool includeSurfaceArea, hasInitializedKernels;
535
    int computeBornSumThreads, gkForceThreads, chainRuleThreads, ediffThreads;
536
    AmoebaMultipoleForce::PolarizationType polarizationType;
537
    std::map<std::string, std::string> defines;
538
539
540
541
542
543
544
545
546
    CudaArray params;
    CudaArray bornSum;
    CudaArray bornRadii;
    CudaArray bornForce;
    CudaArray field;
    CudaArray inducedField;
    CudaArray inducedFieldPolar;
    CudaArray inducedDipoleS;
    CudaArray inducedDipolePolarS;
547
    CUfunction computeBornSumKernel, reduceBornSumKernel, surfaceAreaKernel, gkForceKernel, chainRuleKernel, ediffKernel;
548
549
550
551
552
553
554
};

/**
 * This kernel is invoked to calculate the vdw forces acting on the system and the energy of the system.
 */
class CudaCalcAmoebaVdwForceKernel : public CalcAmoebaVdwForceKernel {
public:
555
    CudaCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system);
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
    ~CudaCalcAmoebaVdwForceKernel();
    /**
     * 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 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);
573
574
575
576
577
578
579
    /**
     * 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);
580
581
582
private:
    class ForceInfo;
    CudaContext& cu;
583
    const System& system;
584
    bool hasInitializedNonbonded;
585
    double dispersionCoefficient;
586
587
588
589
590
    CudaArray sigmaEpsilon;
    CudaArray bondReductionAtoms;
    CudaArray bondReductionFactors;
    CudaArray tempPosq;
    CudaArray tempForces;
591
592
    CudaNonbondedUtilities* nonbonded;
    CUfunction prepareKernel, spreadKernel;
593
594
595
596
597
598
599
};

/**
 * This kernel is invoked to calculate the WCA dispersion forces acting on the system and the energy of the system.
 */
class CudaCalcAmoebaWcaDispersionForceKernel : public CalcAmoebaWcaDispersionForceKernel {
public:
600
    CudaCalcAmoebaWcaDispersionForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system);
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
    /**
     * 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);
617
618
619
620
621
622
623
    /**
     * 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);
624
625
626
private:
    class ForceInfo;
    CudaContext& cu;
627
    const System& system;
628
    double totalMaximumDispersionEnergy;
629
    CudaArray radiusEpsilon;
630
    CUfunction forceKernel;
631
632
633
634
635
};

} // namespace OpenMM

#endif /*AMOEBA_OPENMM_CUDAKERNELS_H*/