CudaNonbondedUtilities.h 18.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#ifndef OPENMM_CUDANONBONDEDUTILITIES_H_
#define OPENMM_CUDANONBONDEDUTILITIES_H_

/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
12
 * Portions copyright (c) 2009-2022 Stanford University and the Authors.      *
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
 * Authors: 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/System.h"
31
#include "CudaArray.h"
32
#include "CudaExpressionUtilities.h"
33
34
#include "openmm/common/NonbondedUtilities.h"
#include <cuda.h>
35
36
37
38
39
#include <sstream>
#include <string>
#include <vector>

namespace OpenMM {
40
    
41
class CudaContext;
42
class CudaSort;
43
44
45

/**
 * This class provides a generic interface for calculating nonbonded interactions.  It does this in two
46
 * ways.  First, it can be used to create kernels that evaluate nonbonded interactions.  Clients
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
 * only need to provide the code for evaluating a single interaction and the list of parameters it depends on.
 * A complete kernel is then synthesized using an appropriate algorithm to evaluate all interactions on all
 * atoms.
 *
 * Second, this class itself creates and invokes a single "default" interaction kernel, allowing several
 * different forces to be evaluated at once for greater efficiency.  Call addInteraction() and addParameter()
 * to add interactions to this default kernel.
 *
 * During each force or energy evaluation, the following sequence of steps takes place:
 *
 * 1. Data structures (e.g. neighbor lists) are calculated to allow nonbonded interactions to be evaluated
 * quickly.
 *
 * 2. calcForcesAndEnergy() is called on each ForceImpl in the System.
 *
 * 3. Finally, the default interaction kernel is invoked to calculate all interactions that were added
 * to it.
 *
 * This sequence means that the default interaction kernel may depend on quantities that were calculated
 * by ForceImpls during calcForcesAndEnergy().
 */

69
class OPENMM_EXPORT_COMMON CudaNonbondedUtilities : public NonbondedUtilities  {
70
71
public:
    class ParameterInfo;
72
73
    CudaNonbondedUtilities(CudaContext& context);
    ~CudaNonbondedUtilities();
74
75
76
77
78
79
80
81
82
83
    /**
     * Add a nonbonded interaction to be evaluated by the default interaction kernel.
     *
     * @param usesCutoff       specifies whether a cutoff should be applied to this interaction
     * @param usesPeriodic     specifies whether periodic boundary conditions should be applied to this interaction
     * @param usesExclusions   specifies whether this interaction uses exclusions.  If this is true, it must have identical exclusions to every other interaction.
     * @param cutoffDistance   the cutoff distance for this interaction (ignored if usesCutoff is false)
     * @param exclusionList    for each atom, specifies the list of other atoms whose interactions should be excluded
     * @param kernel           the code to evaluate the interaction
     * @param forceGroup       the force group in which the interaction should be calculated
84
85
     * @param useNeighborList  specifies whether a neighbor list should be used to optimize this interaction.  This should
     *                         be viewed as only a suggestion.  Even when it is false, a neighbor list may be used anyway.
86
     */
87
    void addInteraction(bool usesCutoff, bool usesPeriodic, bool usesExclusions, double cutoffDistance, const std::vector<std::vector<int> >& exclusionList, const std::string& kernel, int forceGroup, bool useNeighborList=true);
88
89
90
    /**
     * Add a nonbonded interaction to be evaluated by the default interaction kernel.
     *
91
92
93
94
95
96
97
     * @param usesCutoff       specifies whether a cutoff should be applied to this interaction
     * @param usesPeriodic     specifies whether periodic boundary conditions should be applied to this interaction
     * @param usesExclusions   specifies whether this interaction uses exclusions.  If this is true, it must have identical exclusions to every other interaction.
     * @param cutoffDistance   the cutoff distance for this interaction (ignored if usesCutoff is false)
     * @param exclusionList    for each atom, specifies the list of other atoms whose interactions should be excluded
     * @param kernel           the code to evaluate the interaction
     * @param forceGroup       the force group in which the interaction should be calculated
98
99
     * @param useNeighborList  specifies whether a neighbor list should be used to optimize this interaction.  This should
     *                         be viewed as only a suggestion.  Even when it is false, a neighbor list may be used anyway.
100
     * @param supportsPairList specifies whether this interaction can work with a neighbor list that uses a separate pair list
101
     */
102
    void addInteraction(bool usesCutoff, bool usesPeriodic, bool usesExclusions, double cutoffDistance, const std::vector<std::vector<int> >& exclusionList, const std::string& kernel, int forceGroup, bool useNeighborList, bool supportsPairList);
103
104
105
    /**
     * Add a per-atom parameter that the default interaction kernel may depend on.
     */
106
107
108
109
110
111
    void addParameter(ComputeParameterInfo parameter);
    /**
     * Add a per-atom parameter that the default interaction kernel may depend on.
     * 
     * @deprecated Use the version that takes a ComputeParameterInfo instead.
     */
112
113
114
115
    void addParameter(const ParameterInfo& parameter);
    /**
     * Add an array (other than a per-atom parameter) that should be passed as an argument to the default interaction kernel.
     */
116
117
118
119
120
121
    void addArgument(ComputeParameterInfo parameter);
    /**
     * Add an array (other than a per-atom parameter) that should be passed as an argument to the default interaction kernel.
     * 
     * @deprecated Use the version that takes a ComputeParameterInfo instead.
     */
122
    void addArgument(const ParameterInfo& parameter);
123
124
125
126
127
128
129
130
131
    /**
     * Register that the interaction kernel will be computing the derivative of the potential energy
     * with respect to a parameter.
     * 
     * @param param   the name of the parameter
     * @return the variable that will be used to accumulate the derivative.  Any code you pass to addInteraction() should
     * add its contributions to this variable.
     */
    std::string addEnergyParameterDerivative(const std::string& param);
132
133
134
135
136
137
138
139
140
141
    /**
     * Specify the list of exclusions that an interaction outside the default kernel will depend on.
     * 
     * @param exclusionList  for each atom, specifies the list of other atoms whose interactions should be excluded
     */
    void requestExclusions(const std::vector<std::vector<int> >& exclusionList);
    /**
     * Initialize this object in preparation for a simulation.
     */
    void initialize(const System& system);
142
143
144
145
146
147
    /**
     * Get the number of force buffers required for nonbonded forces.
     */
    int getNumForceBuffers() const {
        return 0;
    }
148
149
150
151
152
153
154
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 number of energy buffers required for nonbonded forces.
     */
    int getNumEnergyBuffers() {
        return numForceThreadBlocks*forceThreadBlockSize;
    }
    /**
     * Get whether a cutoff is being used.
     */
    bool getUseCutoff() {
        return useCutoff;
    }
    /**
     * Get whether periodic boundary conditions are being used.
     */
    bool getUsePeriodic() {
        return usePeriodic;
    }
    /**
     * Get the number of work groups used for computing nonbonded forces.
     */
    int getNumForceThreadBlocks() {
        return numForceThreadBlocks;
    }
    /**
     * Get the size of each work group used for computing nonbonded forces.
     */
    int getForceThreadBlockSize() {
        return forceThreadBlockSize;
    }
    /**
179
     * Get the maximum cutoff distance used by any force group.
180
     */
181
    double getMaxCutoffDistance();
182
183
184
185
186
    /**
     * Given a nonbonded cutoff, get the padded cutoff distance used in computing
     * the neighbor list.
     */
    double padCutoff(double cutoff);
187
188
189
    /**
     * Prepare to compute interactions.  This updates the neighbor list.
     */
190
    void prepareInteractions(int forceGroups);
191
192
    /**
     * Compute the nonbonded interactions.
193
194
195
196
     * 
     * @param forceGroups    the flags specifying which force groups to include
     * @param includeForces  whether to compute forces
     * @param includeEnergy  whether to compute the potential energy
197
     */
198
    void computeInteractions(int forceGroups, bool includeForces, bool includeEnergy);
199
200
    /**
     * Check to see if the neighbor list arrays are large enough, and make them bigger if necessary.
201
202
     *
     * @return true if the neighbor list needed to be enlarged.
203
     */
204
    bool updateNeighborListSize();
205
206
207
208
    /**
     * Get the array containing the center of each atom block.
     */
    CudaArray& getBlockCenters() {
Peter Eastman's avatar
Peter Eastman committed
209
        return blockCenter;
210
211
212
213
214
    }
    /**
     * Get the array containing the dimensions of each atom block.
     */
    CudaArray& getBlockBoundingBoxes() {
Peter Eastman's avatar
Peter Eastman committed
215
        return blockBoundingBox;
216
217
218
219
220
    }
    /**
     * Get the array whose first element contains the number of tiles with interactions.
     */
    CudaArray& getInteractionCount() {
Peter Eastman's avatar
Peter Eastman committed
221
        return interactionCount;
222
223
224
225
226
    }
    /**
     * Get the array containing tiles with interactions.
     */
    CudaArray& getInteractingTiles() {
Peter Eastman's avatar
Peter Eastman committed
227
        return interactingTiles;
228
229
    }
    /**
230
     * Get the array containing the atoms in each tile with interactions.
231
     */
232
    CudaArray& getInteractingAtoms() {
Peter Eastman's avatar
Peter Eastman committed
233
        return interactingAtoms;
234
    }
235
236
237
238
    /**
     * Get the array containing single pairs in the neighbor list.
     */
    CudaArray& getSinglePairs() {
Peter Eastman's avatar
Peter Eastman committed
239
        return singlePairs;
240
    }
241
242
243
244
    /**
     * Get the array containing exclusion flags.
     */
    CudaArray& getExclusions() {
Peter Eastman's avatar
Peter Eastman committed
245
        return exclusions;
246
    }
247
248
249
250
    /**
     * Get the array containing tiles with exclusions.
     */
    CudaArray& getExclusionTiles() {
Peter Eastman's avatar
Peter Eastman committed
251
        return exclusionTiles;
252
    }
253
254
255
256
    /**
     * Get the array containing the index into the exclusion array for each tile.
     */
    CudaArray& getExclusionIndices() {
Peter Eastman's avatar
Peter Eastman committed
257
        return exclusionIndices;
258
259
260
261
262
    }
    /**
     * Get the array listing where the exclusion data starts for each row.
     */
    CudaArray& getExclusionRowIndices() {
Peter Eastman's avatar
Peter Eastman committed
263
        return exclusionRowIndices;
264
    }
265
266
267
268
269
270
271
    /**
     * Get the array containing a flag for whether the neighbor list was rebuilt
     * on the most recent call to prepareInteractions().
     */
    CudaArray& getRebuildNeighborList() {
        return rebuildNeighborList;
    }
272
273
274
275
276
277
278
279
280
281
282
283
284
    /**
     * Get the index of the first tile this context is responsible for processing.
     */
    int getStartTileIndex() const {
        return startTileIndex;
    }
    /**
     * Get the total number of tiles this context is responsible for processing.
     */
    int getNumTiles() const {
        return numTiles;
    }
    /**
285
286
287
288
289
290
291
292
293
     * Set whether to add padding to the cutoff distance when building the neighbor list.
     * This increases the size of the neighbor list (and thus the cost of computing interactions),
     * but also means we don't need to rebuild it every time step.  The default value is true,
     * since usually this improves performance.  For very expensive interactions, however,
     * it may be better to set this to false.
     */
    void setUsePadding(bool padding);
    /**
     * Set the range of atom blocks and tiles that should be processed by this context.
294
     */
295
    void setAtomBlockRange(double startFraction, double endFraction);
296
297
298
299
300
301
302
303
304
305
    /**
     * Create a Kernel for evaluating a nonbonded interaction.  Cutoffs and periodic boundary conditions
     * are assumed to be the same as those for the default interaction Kernel, since this kernel will use
     * the same neighbor list.
     * 
     * @param source        the source code for evaluating the force and energy
     * @param params        the per-atom parameters this kernel may depend on
     * @param arguments     arrays (other than per-atom parameters) that should be passed as arguments to the kernel
     * @param useExclusions specifies whether exclusions are applied to this interaction
     * @param isSymmetric   specifies whether the interaction is symmetric
306
     * @param groups        the set of force groups this kernel is for
307
308
     * @param includeForces whether this kernel should compute forces
     * @param includeEnergy whether this kernel should compute potential energy
309
     */
310
    CUfunction createInteractionKernel(const std::string& source, std::vector<ParameterInfo>& params, std::vector<ParameterInfo>& arguments, bool useExclusions, bool isSymmetric, int groups, bool includeForces, bool includeEnergy);
311
312
313
314
315
316
    /**
     * Create the set of kernels that will be needed for a particular combination of force groups.
     * 
     * @param groups    the set of force groups
     */
    void createKernelsForGroups(int groups);
peastman's avatar
peastman committed
317
318
319
320
321
    /**
     * Set the source code for the main kernel.  This defaults to the content of nonbonded.cu.  It only needs to be
     * changed in very unusual circumstances.
     */
    void setKernelSource(const std::string& source);
322
private:
323
    class KernelSet;
324
    class BlockSortTrait;
325
    CudaContext& context;
326
    std::map<int, KernelSet> groupKernels;
Peter Eastman's avatar
Peter Eastman committed
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
    CudaArray exclusionTiles;
    CudaArray exclusions;
    CudaArray exclusionIndices;
    CudaArray exclusionRowIndices;
    CudaArray interactingTiles;
    CudaArray interactingAtoms;
    CudaArray interactionCount;
    CudaArray singlePairs;
    CudaArray singlePairCount;
    CudaArray blockCenter;
    CudaArray blockBoundingBox;
    CudaArray sortedBlocks;
    CudaArray sortedBlockCenter;
    CudaArray sortedBlockBoundingBox;
    CudaArray oldPositions;
    CudaArray rebuildNeighborList;
343
    CudaSort* blockSorter;
344
    CUevent downloadCountEvent;
345
    unsigned int* pinnedCountBuffer;
346
    std::vector<void*> forceArgs, findBlockBoundsArgs, sortBoxDataArgs, findInteractingBlocksArgs;
347
348
349
    std::vector<std::vector<int> > atomExclusions;
    std::vector<ParameterInfo> parameters;
    std::vector<ParameterInfo> arguments;
350
    std::vector<std::string> energyParameterDerivatives;
351
352
353
    std::map<int, double> groupCutoff;
    std::map<int, std::string> groupKernelSource;
    double lastCutoff;
354
    bool useCutoff, usePeriodic, anyExclusions, usePadding, useNeighborList, forceRebuildNeighborList, canUsePairList;
355
    int startTileIndex, startBlockIndex, numBlocks, maxExclusions, numForceThreadBlocks, forceThreadBlockSize, numAtoms, groupFlags;
356
    unsigned int maxTiles, maxSinglePairs, tilesAfterReorder;
357
    long long numTiles;
peastman's avatar
peastman committed
358
    std::string kernelSource;
359
360
361
362
363
364
365
366
367
368
};

/**
 * This class stores the kernels to execute for a set of force groups.
 */

class CudaNonbondedUtilities::KernelSet {
public:
    bool hasForces;
    double cutoffDistance;
369
370
    std::string source;
    CUfunction forceKernel, energyKernel, forceEnergyKernel;
371
372
373
374
    CUfunction findBlockBoundsKernel;
    CUfunction sortBoxDataKernel;
    CUfunction findInteractingBlocksKernel;
    CUfunction findInteractionsWithinBlocksKernel;
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
};

/**
 * This class stores information about a per-atom parameter that may be used in a nonbonded kernel.
 */

class CudaNonbondedUtilities::ParameterInfo {
public:
    /**
     * Create a ParameterInfo object.
     *
     * @param name           the name of the parameter
     * @param type           the data type of the parameter's components
     * @param numComponents  the number of components in the parameter
     * @param size           the size of the parameter in bytes
     * @param memory         the memory containing the parameter values
peastman's avatar
peastman committed
391
     * @param constant       whether the memory should be marked as constant
392
     */
peastman's avatar
peastman committed
393
394
    ParameterInfo(const std::string& name, const std::string& componentType, int numComponents, int size, CUdeviceptr memory, bool constant=true) :
            name(name), componentType(componentType), numComponents(numComponents), size(size), memory(memory), constant(constant) {
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
        if (numComponents == 1)
            type = componentType;
        else {
            std::stringstream s;
            s << componentType << numComponents;
            type = s.str();
        }
    }
    const std::string& getName() const {
        return name;
    }
    const std::string& getComponentType() const {
        return componentType;
    }
    const std::string& getType() const {
        return type;
    }
    int getNumComponents() const {
        return numComponents;
    }
    int getSize() const {
        return size;
    }
418
    CUdeviceptr& getMemory() {
419
420
        return memory;
    }
peastman's avatar
peastman committed
421
422
423
    bool isConstant() const {
        return constant;
    }
424
425
426
427
428
429
private:
    std::string name;
    std::string componentType;
    std::string type;
    int size, numComponents;
    CUdeviceptr memory;
peastman's avatar
peastman committed
430
    bool constant;
431
432
433
434
435
};

} // namespace OpenMM

#endif /*OPENMM_CUDANONBONDEDUTILITIES_H_*/