CudaNonbondedUtilities.h 14.8 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.               *
 *                                                                            *
Peter Eastman's avatar
Peter Eastman committed
12
 * Portions copyright (c) 2009-2025 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
#include "openmm/common/ComputeSort.h"
34
35
#include "openmm/common/NonbondedUtilities.h"
#include <cuda.h>
36
37
38
39
40
#include <sstream>
#include <string>
#include <vector>

namespace OpenMM {
41
    
42
class CudaContext;
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
public:
71
72
73
74
75
    CudaNonbondedUtilities(CudaContext& context);
    ~CudaNonbondedUtilities();
    /**
     * Add a nonbonded interaction to be evaluated by the default interaction kernel.
     *
76
77
78
79
80
81
82
     * @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
83
84
     * @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.
85
     * @param supportsPairList specifies whether this interaction can work with a neighbor list that uses a separate pair list
86
     */
87
88
89
    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, bool supportsPairList=false);
90
91
92
    /**
     * Add a per-atom parameter that the default interaction kernel may depend on.
     */
93
    void addParameter(ComputeParameterInfo parameter);
94
95
96
    /**
     * Add an array (other than a per-atom parameter) that should be passed as an argument to the default interaction kernel.
     */
97
    void addArgument(ComputeParameterInfo parameter);
98
99
100
101
102
103
104
105
106
    /**
     * 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);
107
108
109
110
111
112
113
114
115
116
    /**
     * 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);
117
118
119
120
121
122
    /**
     * Get the number of force buffers required for nonbonded forces.
     */
    int getNumForceBuffers() const {
        return 0;
    }
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
    /**
     * 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;
    }
    /**
154
     * Get the maximum cutoff distance used by any force group.
155
     */
156
    double getMaxCutoffDistance();
157
158
159
160
161
    /**
     * Given a nonbonded cutoff, get the padded cutoff distance used in computing
     * the neighbor list.
     */
    double padCutoff(double cutoff);
162
163
164
    /**
     * Prepare to compute interactions.  This updates the neighbor list.
     */
165
    void prepareInteractions(int forceGroups);
166
167
    /**
     * Compute the nonbonded interactions.
168
169
170
171
     * 
     * @param forceGroups    the flags specifying which force groups to include
     * @param includeForces  whether to compute forces
     * @param includeEnergy  whether to compute the potential energy
172
     */
173
    void computeInteractions(int forceGroups, bool includeForces, bool includeEnergy);
174
175
    /**
     * Check to see if the neighbor list arrays are large enough, and make them bigger if necessary.
176
177
     *
     * @return true if the neighbor list needed to be enlarged.
178
     */
179
    bool updateNeighborListSize();
180
181
182
183
    /**
     * Get the array containing the center of each atom block.
     */
    CudaArray& getBlockCenters() {
Peter Eastman's avatar
Peter Eastman committed
184
        return blockCenter;
185
186
187
188
189
    }
    /**
     * Get the array containing the dimensions of each atom block.
     */
    CudaArray& getBlockBoundingBoxes() {
Peter Eastman's avatar
Peter Eastman committed
190
        return blockBoundingBox;
191
192
193
194
195
    }
    /**
     * Get the array whose first element contains the number of tiles with interactions.
     */
    CudaArray& getInteractionCount() {
Peter Eastman's avatar
Peter Eastman committed
196
        return interactionCount;
197
198
199
200
201
    }
    /**
     * Get the array containing tiles with interactions.
     */
    CudaArray& getInteractingTiles() {
Peter Eastman's avatar
Peter Eastman committed
202
        return interactingTiles;
203
204
    }
    /**
205
     * Get the array containing the atoms in each tile with interactions.
206
     */
207
    CudaArray& getInteractingAtoms() {
Peter Eastman's avatar
Peter Eastman committed
208
        return interactingAtoms;
209
    }
210
211
212
213
    /**
     * Get the array containing single pairs in the neighbor list.
     */
    CudaArray& getSinglePairs() {
Peter Eastman's avatar
Peter Eastman committed
214
        return singlePairs;
215
    }
216
217
218
219
    /**
     * Get the array containing exclusion flags.
     */
    CudaArray& getExclusions() {
Peter Eastman's avatar
Peter Eastman committed
220
        return exclusions;
221
    }
222
223
224
225
    /**
     * Get the array containing tiles with exclusions.
     */
    CudaArray& getExclusionTiles() {
Peter Eastman's avatar
Peter Eastman committed
226
        return exclusionTiles;
227
    }
228
229
230
231
    /**
     * Get the array containing the index into the exclusion array for each tile.
     */
    CudaArray& getExclusionIndices() {
Peter Eastman's avatar
Peter Eastman committed
232
        return exclusionIndices;
233
234
235
236
237
    }
    /**
     * Get the array listing where the exclusion data starts for each row.
     */
    CudaArray& getExclusionRowIndices() {
Peter Eastman's avatar
Peter Eastman committed
238
        return exclusionRowIndices;
239
    }
240
241
242
243
244
245
246
    /**
     * Get the array containing a flag for whether the neighbor list was rebuilt
     * on the most recent call to prepareInteractions().
     */
    CudaArray& getRebuildNeighborList() {
        return rebuildNeighborList;
    }
247
248
249
250
251
252
253
254
255
256
257
258
259
    /**
     * 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;
    }
    /**
260
261
262
263
264
265
266
267
268
     * 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.
269
     */
270
    void setAtomBlockRange(double startFraction, double endFraction);
271
272
273
274
275
276
277
278
279
280
    /**
     * 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
281
     * @param groups        the set of force groups this kernel is for
282
283
     * @param includeForces whether this kernel should compute forces
     * @param includeEnergy whether this kernel should compute potential energy
284
     */
285
    CUfunction createInteractionKernel(const std::string& source, std::vector<ComputeParameterInfo>& params, std::vector<ComputeParameterInfo>& arguments, bool useExclusions, bool isSymmetric, int groups, bool includeForces, bool includeEnergy);
286
287
288
289
290
291
    /**
     * 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
292
293
294
295
296
    /**
     * 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);
297
private:
298
    class KernelSet;
299
    class BlockSortTrait;
300
    void initParamArgs();
301
    CudaContext& context;
302
    std::map<int, KernelSet> groupKernels;
Peter Eastman's avatar
Peter Eastman committed
303
304
305
306
307
308
309
310
311
312
313
314
315
316
    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;
317
    CudaArray blockSizeRange;
318
319
    CudaArray largeBlockCenter;
    CudaArray largeBlockBoundingBox;
Peter Eastman's avatar
Peter Eastman committed
320
321
    CudaArray oldPositions;
    CudaArray rebuildNeighborList;
322
    ComputeSort blockSorter;
323
    CUevent downloadCountEvent;
324
    unsigned int* pinnedCountBuffer;
325
    std::vector<void*> forceArgs, findBlockBoundsArgs, computeSortKeysArgs, sortBoxDataArgs, findInteractingBlocksArgs;
326
    std::vector<std::vector<int> > atomExclusions;
327
328
    std::vector<ComputeParameterInfo> parameters;
    std::vector<ComputeParameterInfo> arguments;
329
    std::vector<std::string> energyParameterDerivatives;
330
331
    std::map<int, double> groupCutoff;
    std::map<int, std::string> groupKernelSource;
Peter Eastman's avatar
Peter Eastman committed
332
    double maxCutoff;
333
334
    bool useCutoff, usePeriodic, anyExclusions, usePadding, useNeighborList, forceRebuildNeighborList, canUsePairList, useLargeBlocks, hasInitializedParams;
    int startTileIndex, startBlockIndex, numBlocks, maxExclusions, numForceThreadBlocks, forceThreadBlockSize, numAtoms, groupFlags, numBlockSizes, paramStartIndex;
335
    unsigned int maxTiles, maxSinglePairs, tilesAfterReorder;
336
    long long numTiles;
peastman's avatar
peastman committed
337
    std::string kernelSource;
338
339
340
341
342
343
344
345
346
};

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

class CudaNonbondedUtilities::KernelSet {
public:
    bool hasForces;
347
348
    std::string source;
    CUfunction forceKernel, energyKernel, forceEnergyKernel;
349
    CUfunction findBlockBoundsKernel;
350
    CUfunction computeSortKeysKernel;
351
352
353
    CUfunction sortBoxDataKernel;
    CUfunction findInteractingBlocksKernel;
    CUfunction findInteractionsWithinBlocksKernel;
354
355
356
357
358
};

} // namespace OpenMM

#endif /*OPENMM_CUDANONBONDEDUTILITIES_H_*/