OpenCLNonbondedUtilities.cpp 18.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
/* -------------------------------------------------------------------------- *
 *                                   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.               *
 *                                                                            *
 * Portions copyright (c) 2009 Stanford University and the Authors.           *
 * 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 "OpenCLNonbondedUtilities.h"
#include "OpenCLArray.h"
29
#include "OpenCLCompact.h"
30
#include "OpenCLKernelSources.h"
31
32
33
34
35
36
#include <map>

using namespace OpenMM;
using namespace std;

OpenCLNonbondedUtilities::OpenCLNonbondedUtilities(OpenCLContext& context) : context(context), cutoff(-1.0), useCutoff(false),
37
38
39
40
41
42
43
44
45
46
47
48
        numForceBuffers(0), tiles(NULL), exclusionIndex(NULL), exclusions(NULL), interactingTiles(NULL), interactionFlags(NULL),
        interactionCount(NULL), blockCenter(NULL), blockBoundingBox(NULL), compact(NULL) {
    // Decide how many force buffers to use.

    forceBufferPerAtomBlock = false;
    numForceBuffers = context.getNumThreadBlocks()*OpenCLContext::ThreadBlockSize/OpenCLContext::TileSize;
    if (numForceBuffers >= context.getNumAtomBlocks()) {
        // For small systems, it is more efficient to have one force buffer per block of 32 atoms instead of one per warp.

        forceBufferPerAtomBlock = true;
        numForceBuffers = context.getNumAtomBlocks();
    }
49
50
51
52
53
54
55
56
57
}

OpenCLNonbondedUtilities::~OpenCLNonbondedUtilities() {
    if (tiles != NULL)
        delete tiles;
    if (exclusionIndex != NULL)
        delete exclusionIndex;
    if (exclusions != NULL)
        delete exclusions;
58
59
60
61
62
63
64
65
66
67
68
69
    if (interactingTiles != NULL)
        delete interactingTiles;
    if (interactionFlags != NULL)
        delete interactionFlags;
    if (interactionCount != NULL)
        delete interactionCount;
    if (blockCenter != NULL)
        delete blockCenter;
    if (blockBoundingBox != NULL)
        delete blockBoundingBox;
    if (compact != NULL)
        delete compact;
70
71
}

72
void OpenCLNonbondedUtilities::addInteraction(bool usesCutoff, bool usesPeriodic, bool usesExclusions, double cutoffDistance, const vector<vector<int> >& exclusionList, const string& kernel) {
73
74
75
76
77
78
79
    if (cutoff != -1.0) {
        if (usesCutoff != useCutoff)
            throw OpenMMException("All Forces must agree on whether to use a cutoff");
        if (usesPeriodic != usePeriodic)
            throw OpenMMException("All Forces must agree on whether to use periodic boundary conditions");
        if (cutoffDistance != cutoff)
            throw OpenMMException("All Forces must use the same cutoff distance");
80
81
    }
    if (usesExclusions && atomExclusions.size() != 0) {
82
        bool sameExclusions = (exclusionList.size() == atomExclusions.size());
83
        for (int i = 0; i < (int) exclusionList.size() && sameExclusions; i++) {
84
85
            if (exclusionList[i].size() != atomExclusions[i].size())
                sameExclusions = false;
86
            for (int j = 0; j < (int) exclusionList[i].size(); j++)
87
88
89
90
91
92
                if (exclusionList[i][j] != atomExclusions[i][j])
                    sameExclusions = false;
        }
        if (!sameExclusions)
            throw OpenMMException("All Forces must have identical exceptions");
    }
93
94
95
96
97
98
    useCutoff = usesCutoff;
    usePeriodic = usesPeriodic;
    cutoff = cutoffDistance;
    kernelSource += kernel+"\n";
    if (usesExclusions)
        atomExclusions = exclusionList;
99
100
}

101
102
void OpenCLNonbondedUtilities::addParameter(const ParameterInfo& parameter) {
    parameters.push_back(parameter);
103
104
}

105
106
107
108
void OpenCLNonbondedUtilities::addArgument(const ParameterInfo& parameter) {
    arguments.push_back(parameter);
}

109
110
111
112
void OpenCLNonbondedUtilities::initialize(const System& system) {
    if (cutoff == -1.0)
        return; // There are no nonbonded interactions in the System.
    
113
114
115
116
    if (atomExclusions.size() == 0) {
        // No exclusions were specifically requested, so just mark every atom as not interacting with itself.
        
        atomExclusions.resize(context.getNumAtoms());
117
        for (int i = 0; i < (int) atomExclusions.size(); i++)
118
119
120
            atomExclusions[i].push_back(i);
    }

121
122
123
124
125
126
127
    // Create the list of tiles.

    int numAtomBlocks = context.getNumAtomBlocks();
    int numTiles = numAtomBlocks*(numAtomBlocks+1)/2;
    tiles = new OpenCLArray<cl_uint>(context, numTiles, "tiles");
    vector<cl_uint> tileVec(tiles->getSize());
    unsigned int count = 0;
128
129
    for (unsigned int y = 0; y < (unsigned int) numAtomBlocks; y++)
        for (unsigned int x = y; x < (unsigned int) numAtomBlocks; x++)
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
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
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
            tileVec[count++] = (x << 17) | (y << 2);

    // Mark which tiles have exclusions.

    for (int atom1 = 0; atom1 < (int) atomExclusions.size(); ++atom1) {
        int x = atom1/OpenCLContext::TileSize;
        for (int j = 0; j < (int) atomExclusions[atom1].size(); ++j) {
            int atom2 = atomExclusions[atom1][j];
            int y = atom2/OpenCLContext::TileSize;
            int index = (x > y ? x+y*numAtomBlocks-y*(y+1)/2 : y+x*numAtomBlocks-x*(x+1)/2);
            tileVec[index] |= 1;
        }
    }
    if (context.getPaddedNumAtoms() > context.getNumAtoms()) {
        int lastTile = context.getNumAtoms()/OpenCLContext::TileSize;
        for (int i = 0; i < numTiles; ++i) {
            int x = tileVec[i]>>17;
            int y = (tileVec[i]>>2)&0x7FFF;
            if (x == lastTile || y == lastTile)
                tileVec[i] |= 1;
        }
    }

    // Build a list of indices for the tiles with exclusions.

    exclusionIndex = new OpenCLArray<cl_uint>(context, numTiles, "exclusionIndex");
    vector<cl_uint> exclusionIndexVec(exclusionIndex->getSize());
    int numWithExclusions = 0;
    for (int i = 0; i < numTiles; ++i)
        if ((tileVec[i]&1) == 1)
            exclusionIndexVec[i] = (numWithExclusions++)*OpenCLContext::TileSize;

    // Record the exclusion data.

    exclusions = new OpenCLArray<cl_uint>(context, numWithExclusions*OpenCLContext::TileSize, "exclusions");
    vector<cl_uint> exclusionVec(exclusions->getSize());
    for (int i = 0; i < exclusions->getSize(); ++i)
        exclusionVec[i] = 0xFFFFFFFF;
    for (int atom1 = 0; atom1 < (int) atomExclusions.size(); ++atom1) {
        int x = atom1/OpenCLContext::TileSize;
        int offset1 = atom1-x*OpenCLContext::TileSize;
        for (int j = 0; j < (int) atomExclusions[atom1].size(); ++j) {
            int atom2 = atomExclusions[atom1][j];
            int y = atom2/OpenCLContext::TileSize;
            int offset2 = atom2-y*OpenCLContext::TileSize;
            if (x > y) {
                int tile = x+y*numAtomBlocks-y*(y+1)/2;
                exclusionVec[exclusionIndexVec[tile]+offset1] &= 0xFFFFFFFF-(1<<offset2);
            }
            else {
                int tile = y+x*numAtomBlocks-x*(x+1)/2;
                exclusionVec[exclusionIndexVec[tile]+offset2] &= 0xFFFFFFFF-(1<<offset1);
            }
        }
    }

    // Mark all interactions that involve a padding atom as being excluded.

    for (int atom1 = context.getNumAtoms(); atom1 < context.getPaddedNumAtoms(); ++atom1) {
        int x = atom1/OpenCLContext::TileSize;
        int offset1 = atom1-x*OpenCLContext::TileSize;
        for (int atom2 = 0; atom2 < context.getPaddedNumAtoms(); ++atom2) {
            int y = atom2/OpenCLContext::TileSize;
            int offset2 = atom2-y*OpenCLContext::TileSize;
            if (x >= y) {
                int tile = x+y*numAtomBlocks-y*(y+1)/2;
                exclusionVec[exclusionIndexVec[tile]+offset1] &= 0xFFFFFFFF-(1<<offset2);
            }
            if (y >= x) {
                int tile = y+x*numAtomBlocks-x*(x+1)/2;
                exclusionVec[exclusionIndexVec[tile]+offset2] &= 0xFFFFFFFF-(1<<offset1);
            }
        }
    }
    atomExclusions.clear(); // We won't use this again, so free the memory it used
    tiles->upload(tileVec);
    exclusions->upload(exclusionVec);
    exclusionIndex->upload(exclusionIndexVec);
208
209
210

    // Record the periodic box size.

211
212
    Vec3 boxVectors[3];
    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
213
    periodicBoxSize = mm_float4((float) boxVectors[0][0], (float) boxVectors[1][1], (float) boxVectors[2][2], 0.0f);
214
215
216
217
218
219
220
221
222
223
224

    // Create data structures for the neighbor list.

    if (useCutoff) {
        interactingTiles = new OpenCLArray<cl_uint>(context, numTiles, "interactingTiles");
        interactionFlags = new OpenCLArray<cl_uint>(context, numTiles, "interactionFlags");
        interactionCount = new OpenCLArray<cl_uint>(context, 1, "interactionCount");
        blockCenter = new OpenCLArray<mm_float4>(context, numAtomBlocks, "blockCenter");
        blockBoundingBox = new OpenCLArray<mm_float4>(context, numAtomBlocks, "blockBoundingBox");
        compact = new OpenCLCompact(context);
    }
225
226
227

    // Create kernels.

228
    forceKernel = createInteractionKernel(kernelSource, parameters, arguments, true);
229
230
231
232
233
234
    if (useCutoff) {
        map<string, string> defines;
        if (forceBufferPerAtomBlock)
            defines["USE_OUTPUT_BUFFER_PER_BLOCK"] = "1";
        if (usePeriodic)
            defines["USE_PERIODIC"] = "1";
235
        cl::Program interactingBlocksProgram = context.createProgram(OpenCLKernelSources::findInteractingBlocks, defines);
236
237
238
239
240
241
242
243
        findBlockBoundsKernel = cl::Kernel(interactingBlocksProgram, "findBlockBounds");
        findBlockBoundsKernel.setArg<cl_int>(0, context.getNumAtoms());
        findBlockBoundsKernel.setArg<mm_float4>(1, periodicBoxSize);
        findBlockBoundsKernel.setArg<cl::Buffer>(2, context.getPosq().getDeviceBuffer());
        findBlockBoundsKernel.setArg<cl::Buffer>(3, blockCenter->getDeviceBuffer());
        findBlockBoundsKernel.setArg<cl::Buffer>(4, blockBoundingBox->getDeviceBuffer());
        findInteractingBlocksKernel = cl::Kernel(interactingBlocksProgram, "findBlocksWithInteractions");
        findInteractingBlocksKernel.setArg<cl_int>(0, tiles->getSize());
244
        findInteractingBlocksKernel.setArg<cl_float>(1, (cl_float) (cutoff*cutoff));
245
246
247
248
249
250
        findInteractingBlocksKernel.setArg<mm_float4>(2, periodicBoxSize);
        findInteractingBlocksKernel.setArg<cl::Buffer>(3, tiles->getDeviceBuffer());
        findInteractingBlocksKernel.setArg<cl::Buffer>(4, blockCenter->getDeviceBuffer());
        findInteractingBlocksKernel.setArg<cl::Buffer>(5, blockBoundingBox->getDeviceBuffer());
        findInteractingBlocksKernel.setArg<cl::Buffer>(6, interactionFlags->getDeviceBuffer());
        findInteractionsWithinBlocksKernel = cl::Kernel(interactingBlocksProgram, "findInteractionsWithinBlocks");
251
        findInteractionsWithinBlocksKernel.setArg<cl_float>(0, (cl_float) (cutoff*cutoff));
252
253
254
255
256
257
258
259
260
        findInteractionsWithinBlocksKernel.setArg<mm_float4>(1, periodicBoxSize);
        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(2, context.getPosq().getDeviceBuffer());
        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(3, interactingTiles->getDeviceBuffer());
        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(4, blockCenter->getDeviceBuffer());
        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(5, blockBoundingBox->getDeviceBuffer());
        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(6, interactionFlags->getDeviceBuffer());
        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(7, interactionCount->getDeviceBuffer());
        findInteractionsWithinBlocksKernel.setArg(8, OpenCLContext::ThreadBlockSize*sizeof(cl_uint), NULL);
    }
261
262
263
264
265
}

void OpenCLNonbondedUtilities::prepareInteractions() {
    if (!useCutoff)
        return;
266
267
268
269
270
271

    // Compute the neighbor list.

    context.executeKernel(findBlockBoundsKernel, context.getNumAtoms());
    context.executeKernel(findInteractingBlocksKernel, context.getNumAtoms());
    compact->compactStream(*interactingTiles, *tiles, *interactionFlags, *interactionCount);
Peter Eastman's avatar
Peter Eastman committed
272
273
    if (context.getSIMDWidth() == 32)
        context.executeKernel(findInteractionsWithinBlocksKernel, context.getNumAtoms());
274
275
276
}

void OpenCLNonbondedUtilities::computeInteractions() {
Peter Eastman's avatar
Bug fix  
Peter Eastman committed
277
278
    if (tiles != NULL)
        context.executeKernel(forceKernel, tiles->getSize()*OpenCLContext::TileSize);
279
280
}

281
cl::Kernel OpenCLNonbondedUtilities::createInteractionKernel(const string& source, const vector<ParameterInfo>& params, const vector<ParameterInfo>& arguments, bool useExclusions) const {
282
283
284
    map<string, string> replacements;
    replacements["COMPUTE_INTERACTION"] = source;
    stringstream args;
285
    for (int i = 0; i < (int) params.size(); i++) {
286
287
288
289
290
291
292
293
294
        args << ", __global ";
        args << params[i].getType();
        args << "* global_";
        args << params[i].getName();
        args << ", __local ";
        args << params[i].getType();
        args << "* local_";
        args << params[i].getName();
    }
295
    for (int i = 0; i < (int) arguments.size(); i++) {
296
297
298
299
300
301
302
303
        if ((arguments[i].getBuffer().getInfo<CL_MEM_FLAGS>() & CL_MEM_READ_ONLY) == 0)
            args << ", __global ";
        else
            args << ", __constant ";
        args << arguments[i].getType();
        args << "* ";
        args << arguments[i].getName();
    }
304
305
    replacements["PARAMETER_ARGUMENTS"] = args.str();
    stringstream loadLocal1;
306
    for (int i = 0; i < (int) params.size(); i++) {
307
308
309
310
        loadLocal1 << "local_";
        loadLocal1 << params[i].getName();
        loadLocal1 << "[get_local_id(0)] = ";
        loadLocal1 << params[i].getName();
311
        loadLocal1 << "1;\n";
312
313
314
    }
    replacements["LOAD_LOCAL_PARAMETERS_FROM_1"] = loadLocal1.str();
    stringstream loadLocal2;
315
    for (int i = 0; i < (int) params.size(); i++) {
316
317
318
319
        loadLocal2 << "local_";
        loadLocal2 << params[i].getName();
        loadLocal2 << "[get_local_id(0)] = global_";
        loadLocal2 << params[i].getName();
320
        loadLocal2 << "[j];\n";
321
322
323
    }
    replacements["LOAD_LOCAL_PARAMETERS_FROM_GLOBAL"] = loadLocal2.str();
    stringstream load1;
324
    for (int i = 0; i < (int) params.size(); i++) {
325
326
327
328
329
        load1 << params[i].getType();
        load1 << " ";
        load1 << params[i].getName();
        load1 << "1 = global_";
        load1 << params[i].getName();
330
        load1 << "[atom1];\n";
331
332
333
    }
    replacements["LOAD_ATOM1_PARAMETERS"] = load1.str();
    stringstream load2j;
334
    for (int i = 0; i < (int) params.size(); i++) {
335
336
337
338
339
        load2j << params[i].getType();
        load2j << " ";
        load2j << params[i].getName();
        load2j << "2 = local_";
        load2j << params[i].getName();
340
        load2j << "[atom2];\n";
341
    }
342
    replacements["LOAD_ATOM2_PARAMETERS"] = load2j.str();
343
344
345
346
347
348
349
350
351
    map<string, string> defines;
    if (forceBufferPerAtomBlock)
        defines["USE_OUTPUT_BUFFER_PER_BLOCK"] = "1";
    if (useCutoff)
        defines["USE_CUTOFF"] = "1";
    if (usePeriodic)
        defines["USE_PERIODIC"] = "1";
    if (useExclusions)
        defines["USE_EXCLUSIONS"] = "1";
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
    stringstream xsize, ysize, zsize, cutoffSquared;
    xsize.precision(8);
    ysize.precision(8);
    zsize.precision(8);
    cutoffSquared.precision(8);
    xsize << scientific << periodicBoxSize.x << "f";
    ysize << scientific << periodicBoxSize.y << "f";
    zsize << scientific << periodicBoxSize.z << "f";
    cutoffSquared << scientific << (cutoff*cutoff) << "f";
    defines["PERIODIC_BOX_SIZE_X"] = xsize.str();
    defines["PERIODIC_BOX_SIZE_Y"] = ysize.str();
    defines["PERIODIC_BOX_SIZE_Z"] = zsize.str();
    defines["CUTOFF_SQUARED"] = cutoffSquared.str();
    stringstream natom, padded;
    natom << context.getNumAtoms();
    padded << context.getPaddedNumAtoms();
    defines["NUM_ATOMS"] = natom.str();
    defines["PADDED_NUM_ATOMS"] = padded.str();
370
371
    string file = (context.getSIMDWidth() == 32 ? OpenCLKernelSources::nonbonded_nvidia : OpenCLKernelSources::nonbonded_default);
    cl::Program program = context.createProgram(context.replaceStrings(file, replacements), defines);
372
373
374
    cl::Kernel kernel(program, "computeNonbonded");

    // Set arguments to the Kernel.
375

376
377
378
379
380
381
382
    kernel.setArg<cl::Buffer>(0, context.getForceBuffers().getDeviceBuffer());
    kernel.setArg<cl::Buffer>(1, context.getEnergyBuffer().getDeviceBuffer());
    kernel.setArg<cl::Buffer>(2, context.getPosq().getDeviceBuffer());
    kernel.setArg<cl::Buffer>(3, exclusions->getDeviceBuffer());
    kernel.setArg<cl::Buffer>(4, exclusionIndex->getDeviceBuffer());
    kernel.setArg(5, OpenCLContext::ThreadBlockSize*sizeof(cl_float4), NULL);
    kernel.setArg(6, OpenCLContext::ThreadBlockSize*sizeof(cl_float4), NULL);
383
384
    kernel.setArg(7, OpenCLContext::ThreadBlockSize*sizeof(cl_float4), NULL);
    int paramBase = 10;
385
    if (useCutoff) {
386
        paramBase = 11;
387
388
389
        kernel.setArg<cl::Buffer>(8, interactingTiles->getDeviceBuffer());
        kernel.setArg<cl::Buffer>(9, interactionFlags->getDeviceBuffer());
        kernel.setArg<cl::Buffer>(10, interactionCount->getDeviceBuffer());
390
    }
391
    else {
392
393
        kernel.setArg<cl::Buffer>(8, tiles->getDeviceBuffer());
        kernel.setArg<cl_uint>(9, tiles->getSize());
394
    }
395
396
397
    for (int i = 0; i < (int) params.size(); i++) {
        kernel.setArg<cl::Buffer>(i*2+paramBase, params[i].getBuffer());
        kernel.setArg(i*2+paramBase+1, OpenCLContext::ThreadBlockSize*params[i].getSize(), NULL);
398
    }
399
400
401
402
    paramBase += 2*params.size();
    for (int i = 0; i < (int) arguments.size(); i++) {
        kernel.setArg<cl::Buffer>(i+paramBase, arguments[i].getBuffer());
    }
403
    return kernel;
404
}