CpuKernels.cpp 69.6 KB
Newer Older
1
2
3
4
5
6
7
8
/* -------------------------------------------------------------------------- *
 *                                   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.               *
 *                                                                            *
9
 * Portions copyright (c) 2013-2025 Stanford University and the Authors.      *
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 * Authors: Peter Eastman                                                     *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

#include "CpuKernels.h"
peastman's avatar
peastman committed
33
#include "ReferenceAngleBondIxn.h"
34
#include "ReferenceBondForce.h"
35
#include "ReferenceConstraints.h"
36
37
#include "ReferenceKernelFactory.h"
#include "ReferenceKernels.h"
38
#include "ReferenceLJCoulomb14.h"
39
40
#include "ReferenceProperDihedralBond.h"
#include "ReferenceRbDihedralBond.h"
41
#include "ReferenceTabulatedFunction.h"
42
#include "openmm/Context.h"
43
#include "openmm/OpenMMException.h"
peastman's avatar
peastman committed
44
#include "openmm/Vec3.h"
45
46
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/NonbondedForceImpl.h"
47
#include "openmm/internal/vectorize.h"
48
#include "lepton/CompiledExpression.h"
49
#include "lepton/CustomFunction.h"
50
#include "lepton/Operation.h"
51
#include "lepton/Parser.h"
52
#include <iostream>
53
#include "lepton/ParsedExpression.h"
54
55
56
57

using namespace OpenMM;
using namespace std;

peastman's avatar
peastman committed
58
static vector<Vec3>& extractPositions(ContextImpl& context) {
59
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
60
    return *data->positions;
61
62
}

peastman's avatar
peastman committed
63
static vector<Vec3>& extractVelocities(ContextImpl& context) {
64
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
65
    return *data->velocities;
66
67
}

peastman's avatar
peastman committed
68
static vector<Vec3>& extractForces(ContextImpl& context) {
69
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
70
    return *data->forces;
71
72
}

peastman's avatar
peastman committed
73
static Vec3& extractBoxSize(ContextImpl& context) {
74
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
75
    return *data->periodicBoxSize;
76
77
}

peastman's avatar
peastman committed
78
static Vec3* extractBoxVectors(ContextImpl& context) {
79
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
80
    return data->periodicBoxVectors;
81
82
}

83
84
static ReferenceConstraints& extractConstraints(ContextImpl& context) {
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
85
    return *data->constraints;
86
87
}

88
89
90
91
92
static const ReferenceVirtualSites& extractVirtualSites(ContextImpl& context) {
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
    return *data->virtualSites;
}

93
94
static map<string, double>& extractEnergyParameterDerivatives(ContextImpl& context) {
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
95
    return *data->energyParameterDerivatives;
96
97
}

98
99
100
101
102
103
104
/**
 * Make sure an expression doesn't use any undefined variables.
 */
static void validateVariables(const Lepton::ExpressionTreeNode& node, const set<string>& variables) {
    const Lepton::Operation& op = node.getOperation();
    if (op.getId() == Lepton::Operation::VARIABLE && variables.find(op.getName()) == variables.end())
        throw OpenMMException("Unknown variable in expression: "+op.getName());
peastman's avatar
peastman committed
105
106
    for (auto& child : node.getChildren())
        validateVariables(child, variables);
107
108
}

109
110
111
112
113
/**
 * Compute the kinetic energy of the system, possibly shifting the velocities in time to account
 * for a leapfrog integrator.
 */
static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& masses, double timeShift) {
peastman's avatar
peastman committed
114
115
116
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
117
118
119
    int numParticles = context.getSystem().getNumParticles();
    
    // Compute the shifted velocities.
120

peastman's avatar
peastman committed
121
    vector<Vec3> shiftedVel(numParticles);
122
123
124
125
126
127
128
129
    for (int i = 0; i < numParticles; ++i) {
        if (masses[i] > 0)
            shiftedVel[i] = velData[i]+forceData[i]*(timeShift/masses[i]);
        else
            shiftedVel[i] = velData[i];
    }
    
    // Apply constraints to them.
130
131
132
133
134
135
136

    if (timeShift != 0) {
        vector<double> inverseMasses(numParticles);
        for (int i = 0; i < numParticles; i++)
            inverseMasses[i] = (masses[i] == 0 ? 0 : 1/masses[i]);
        extractConstraints(context).applyToVelocities(posData, shiftedVel, inverseMasses, 1e-4);
    }
137
138
139
140
141
142
143
144
145
146
    
    // Compute the kinetic energy.
    
    double energy = 0.0;
    for (int i = 0; i < numParticles; ++i)
        if (masses[i] > 0)
            energy += masses[i]*(shiftedVel[i].dot(shiftedVel[i]));
    return 0.5*energy;
}

147
148
149
150
151
152
153
154
155
156
157
158
159
/**
 * Copy particle charges into the fourth element of the posq array.
 */
static void copyChargesToPosq(ContextImpl& context, const vector<float>& charges, int index) {
    CpuPlatform::PlatformData& data = CpuPlatform::getPlatformData(context);
    if (index == data.currentPosqIndex)
        return;
    data.currentPosqIndex = index;
    AlignedArray<float>& posq = data.posq;
    for (int i = 0; i < charges.size(); i++)
        posq[4*i+3] = charges[i];
}

peastman's avatar
peastman committed
160
161
162
163
164
165
166
CpuCalcForcesAndEnergyKernel::CpuCalcForcesAndEnergyKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data, ContextImpl& context) :
        CalcForcesAndEnergyKernel(name, platform), data(data) {
    // Create a Reference platform version of this kernel.
    
    ReferenceKernelFactory referenceFactory;
    referenceKernel = Kernel(referenceFactory.createKernelImpl(name, platform, context));
}
167

peastman's avatar
peastman committed
168
169
170
171
172
173
174
175
176
177
178
179
180
void CpuCalcForcesAndEnergyKernel::initialize(const System& system) {
    referenceKernel.getAs<ReferenceCalcForcesAndEnergyKernel>().initialize(system);
    lastPositions.resize(system.getNumParticles(), Vec3(1e10, 1e10, 1e10));
}

void CpuCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups) {
    referenceKernel.getAs<ReferenceCalcForcesAndEnergyKernel>().beginComputation(context, includeForce, includeEnergy, groups);
    
    // Convert positions to single precision and clear the forces.

    int numParticles = context.getSystem().getNumParticles();
    bool positionsValid = true;
    data.threads.execute([&] (ThreadPool& threads, int threadIndex) {
181
182
183
        // Convert the positions to single precision and apply periodic boundary conditions

        AlignedArray<float>& posq = data.posq;
peastman's avatar
peastman committed
184
185
        vector<Vec3>& posData = extractPositions(context);
        Vec3* boxVectors = extractBoxVectors(context);
186
187
188
        double boxSize[3] = {boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]};
        double invBoxSize[3] = {1/boxVectors[0][0], 1/boxVectors[1][1], 1/boxVectors[2][2]};
        bool triclinic = (boxVectors[0][1] != 0 || boxVectors[0][2] != 0 || boxVectors[1][0] != 0 || boxVectors[1][2] != 0 || boxVectors[2][0] != 0 || boxVectors[2][1] != 0);
189
190
191
192
        int numParticles = context.getSystem().getNumParticles();
        int numThreads = threads.getNumThreads();
        int start = threadIndex*numParticles/numThreads;
        int end = (threadIndex+1)*numParticles/numThreads;
193
194
195
        if (data.isPeriodic) {
            if (triclinic) {
                for (int i = start; i < end; i++) {
peastman's avatar
peastman committed
196
                    Vec3 pos = posData[i];
197
198
199
200
201
202
                    pos -= boxVectors[2]*floor(pos[2]*invBoxSize[2]);
                    pos -= boxVectors[1]*floor(pos[1]*invBoxSize[1]);
                    pos -= boxVectors[0]*floor(pos[0]*invBoxSize[0]);
                    posq[4*i] = (float) pos[0];
                    posq[4*i+1] = (float) pos[1];
                    posq[4*i+2] = (float) pos[2];
203
                }
204
205
206
207
            }
            else {
                for (int i = start; i < end; i++) {
                    for (int j = 0; j < 3; j++) {
peastman's avatar
peastman committed
208
                        double x = posData[i][j];
209
210
211
                        double base = floor(x*invBoxSize[j])*boxSize[j];
                        posq[4*i+j] = (float) (x-base);
                    }
212
                }
213
214
            }
        }
215
216
217
218
219
220
        else
            for (int i = start; i < end; i++) {
                posq[4*i] = (float) posData[i][0];
                posq[4*i+1] = (float) posData[i][1];
                posq[4*i+2] = (float) posData[i][2];
            }
221
222
223
        
        // Check for invalid positions.
        
224
225
        for (int i = 4*start; i < 4*end; i += 4)
            if (posq[i] != posq[i] || posq[i+1] != posq[i+1] || posq[i+2] != posq[i+2])
226
                positionsValid = false;
227
228
229
230
231
232

        // Clear the forces.

        fvec4 zero(0.0f);
        for (int j = 0; j < numParticles; j++)
            zero.store(&data.threadForce[threadIndex][j*4]);
peastman's avatar
peastman committed
233
    });
234
    data.threads.waitForThreads();
peastman's avatar
peastman committed
235
    if (!positionsValid)
236
        throw OpenMMException("Particle coordinate is NaN.  For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#nan");
237
238
239

    // Determine whether we need to recompute the neighbor list.
        
240
    if (data.neighborList != NULL && data.cutoff > 0.0) {
241
242
243
244
245
246
        double padding = data.paddedCutoff-data.cutoff;;
        bool needRecompute = false;
        double closeCutoff2 = 0.25*padding*padding;
        double farCutoff2 = 0.5*padding*padding;
        int maxNumMoved = numParticles/10;
        vector<int> moved;
peastman's avatar
peastman committed
247
        vector<Vec3>& posData = extractPositions(context);
248
        for (int i = 0; i < numParticles; i++) {
peastman's avatar
peastman committed
249
            Vec3 delta = posData[i]-lastPositions[i];
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
            double dist2 = delta.dot(delta);
            if (dist2 > closeCutoff2) {
                moved.push_back(i);
                if (dist2 > farCutoff2 || moved.size() > maxNumMoved) {
                    needRecompute = true;
                    break;
                }
            }
        }
        if (!needRecompute && moved.size() > 0) {
            // Some particles have moved further than half the padding distance.  Look for pairs
            // that are missing from the neighbor list.

            int numMoved = moved.size();
            double cutoff2 = data.cutoff*data.cutoff;
            double paddedCutoff2 = data.paddedCutoff*data.paddedCutoff;
            for (int i = 1; i < numMoved && !needRecompute; i++)
                for (int j = 0; j < i; j++) {
peastman's avatar
peastman committed
268
                    Vec3 delta = posData[moved[i]]-posData[moved[j]];
269
270
271
                    if (delta.dot(delta) < cutoff2) {
                        // These particles should interact.  See if they are in the neighbor list.
                        
peastman's avatar
peastman committed
272
                        Vec3 oldDelta = lastPositions[moved[i]]-lastPositions[moved[j]];
273
274
275
276
277
278
279
280
281
282
283
284
                        if (oldDelta.dot(oldDelta) > paddedCutoff2) {
                            needRecompute = true;
                            break;
                        }
                    }
                }
        }
        if (needRecompute) {
            data.neighborList->computeNeighborList(numParticles, data.posq, data.exclusions, extractBoxVectors(context), data.isPeriodic, data.paddedCutoff, data.threads);
            lastPositions = posData;
        }
    }
285
286
}

287
double CpuCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForce, bool includeEnergy, int groups, bool& valid) {
288
289
    // Sum the forces from all the threads.
    
peastman's avatar
peastman committed
290
291
292
293
294
295
296
    data.threads.execute([&] (ThreadPool& threads, int threadIndex) {
        // Sum the contributions to forces that have been calculated by different threads.
        
        int numParticles = context.getSystem().getNumParticles();
        int numThreads = threads.getNumThreads();
        int start = threadIndex*numParticles/numThreads;
        int end = (threadIndex+1)*numParticles/numThreads;
peastman's avatar
peastman committed
297
        vector<Vec3>& forceData = extractForces(context);
peastman's avatar
peastman committed
298
299
300
301
302
303
304
305
306
        for (int i = start; i < end; i++) {
            fvec4 f(0.0f);
            for (int j = 0; j < numThreads; j++)
                f += fvec4(&data.threadForce[j][4*i]);
            forceData[i][0] += f[0];
            forceData[i][1] += f[1];
            forceData[i][2] += f[2];
        }
    });
307
    data.threads.waitForThreads();
308
    return referenceKernel.getAs<ReferenceCalcForcesAndEnergyKernel>().finishComputation(context, includeForce, includeEnergy, groups, valid);
309
310
}

311
312
313
314
315
316
317
318
319
320
void CpuUpdateStateDataKernel::createCheckpoint(ContextImpl& context, ostream& stream) {
    ReferenceUpdateStateDataKernel::createCheckpoint(context, stream);
    data.random.createCheckpoint(stream);
}

void CpuUpdateStateDataKernel::loadCheckpoint(ContextImpl& context, istream& stream) {
    ReferenceUpdateStateDataKernel::loadCheckpoint(context, stream);
    data.random.loadCheckpoint(stream);
}

peastman's avatar
peastman committed
321
322
void CpuCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
    numAngles = force.getNumAngles();
323
324
    angleIndexArray.resize(numAngles, vector<int>(3));
    angleParamArray.resize(numAngles, vector<double>(2));
peastman's avatar
peastman committed
325
326
327
328
329
330
331
    for (int i = 0; i < numAngles; ++i) {
        int particle1, particle2, particle3;
        double angle, k;
        force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
        angleIndexArray[i][0] = particle1;
        angleIndexArray[i][1] = particle2;
        angleIndexArray[i][2] = particle3;
peastman's avatar
peastman committed
332
333
        angleParamArray[i][0] = angle;
        angleParamArray[i][1] = k;
peastman's avatar
peastman committed
334
335
    }
    bondForce.initialize(system.getNumParticles(), numAngles, 3, angleIndexArray, data.threads);
336
    usePeriodic = force.usesPeriodicBoundaryConditions();
peastman's avatar
peastman committed
337
338
339
}

double CpuCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
340
341
342
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
peastman's avatar
peastman committed
343
    ReferenceAngleBondIxn angleBond;
344
345
    if (usePeriodic)
        angleBond.setPeriodic(extractBoxVectors(context));
peastman's avatar
peastman committed
346
347
348
349
    bondForce.calculateForce(posData, angleParamArray, forceData, includeEnergy ? &energy : NULL, angleBond);
    return energy;
}

350
void CpuCalcHarmonicAngleForceKernel::copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force, int firstAngle, int lastAngle) {
peastman's avatar
peastman committed
351
352
353
354
355
    if (numAngles != force.getNumAngles())
        throw OpenMMException("updateParametersInContext: The number of angles has changed");

    // Record the values.

356
    for (int i = firstAngle; i <= lastAngle; ++i) {
peastman's avatar
peastman committed
357
358
359
360
361
        int particle1, particle2, particle3;
        double angle, k;
        force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
        if (particle1 != angleIndexArray[i][0] || particle2 != angleIndexArray[i][1] || particle3 != angleIndexArray[i][2])
            throw OpenMMException("updateParametersInContext: The set of particles in an angle has changed");
peastman's avatar
peastman committed
362
363
        angleParamArray[i][0] = angle;
        angleParamArray[i][1] = k;
peastman's avatar
peastman committed
364
365
366
    }
}

367
368
void CpuCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
    numTorsions = force.getNumTorsions();
369
370
    torsionIndexArray.resize(numTorsions, vector<int>(4));
    torsionParamArray.resize(numTorsions, vector<double>(3));
371
372
373
374
375
376
377
378
    for (int i = 0; i < numTorsions; ++i) {
        int particle1, particle2, particle3, particle4, periodicity;
        double phase, k;
        force.getTorsionParameters(i, particle1, particle2, particle3, particle4, periodicity, phase, k);
        torsionIndexArray[i][0] = particle1;
        torsionIndexArray[i][1] = particle2;
        torsionIndexArray[i][2] = particle3;
        torsionIndexArray[i][3] = particle4;
peastman's avatar
peastman committed
379
380
381
        torsionParamArray[i][0] = k;
        torsionParamArray[i][1] = phase;
        torsionParamArray[i][2] = periodicity;
382
383
    }
    bondForce.initialize(system.getNumParticles(), numTorsions, 4, torsionIndexArray, data.threads);
384
    usePeriodic = force.usesPeriodicBoundaryConditions();
385
386
387
}

double CpuCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
388
389
390
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
391
    ReferenceProperDihedralBond periodicTorsionBond;
392
393
    if (usePeriodic)
        periodicTorsionBond.setPeriodic(extractBoxVectors(context));
394
395
396
397
    bondForce.calculateForce(posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, periodicTorsionBond);
    return energy;
}

398
void CpuCalcPeriodicTorsionForceKernel::copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& force, int firstTorsion, int lastTorsion) {
399
400
401
402
403
    if (numTorsions != force.getNumTorsions())
        throw OpenMMException("updateParametersInContext: The number of torsions has changed");

    // Record the values.

404
    for (int i = firstTorsion; i <= lastTorsion; ++i) {
405
406
407
408
409
        int particle1, particle2, particle3, particle4, periodicity;
        double phase, k;
        force.getTorsionParameters(i, particle1, particle2, particle3, particle4, periodicity, phase, k);
        if (particle1 != torsionIndexArray[i][0] || particle2 != torsionIndexArray[i][1] || particle3 != torsionIndexArray[i][2] || particle4 != torsionIndexArray[i][3])
            throw OpenMMException("updateParametersInContext: The set of particles in a torsion has changed");
peastman's avatar
peastman committed
410
411
412
        torsionParamArray[i][0] = k;
        torsionParamArray[i][1] = phase;
        torsionParamArray[i][2] = periodicity;
413
414
415
416
417
    }
}

void CpuCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
    numTorsions = force.getNumTorsions();
418
419
    torsionIndexArray.resize(numTorsions, vector<int>(4));
    torsionParamArray.resize(numTorsions, vector<double>(6));
420
421
422
423
424
425
426
427
    for (int i = 0; i < numTorsions; ++i) {
        int particle1, particle2, particle3, particle4;
        double c0, c1, c2, c3, c4, c5;
        force.getTorsionParameters(i, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
        torsionIndexArray[i][0] = particle1;
        torsionIndexArray[i][1] = particle2;
        torsionIndexArray[i][2] = particle3;
        torsionIndexArray[i][3] = particle4;
peastman's avatar
peastman committed
428
429
430
431
432
433
        torsionParamArray[i][0] = c0;
        torsionParamArray[i][1] = c1;
        torsionParamArray[i][2] = c2;
        torsionParamArray[i][3] = c3;
        torsionParamArray[i][4] = c4;
        torsionParamArray[i][5] = c5;
434
435
    }
    bondForce.initialize(system.getNumParticles(), numTorsions, 4, torsionIndexArray, data.threads);
436
    usePeriodic = force.usesPeriodicBoundaryConditions();
437
438
439
}

double CpuCalcRBTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
440
441
442
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
443
    ReferenceRbDihedralBond rbTorsionBond;
444
445
    if (usePeriodic)
        rbTorsionBond.setPeriodic(extractBoxVectors(context));
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
    bondForce.calculateForce(posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, rbTorsionBond);
    return energy;
}

void CpuCalcRBTorsionForceKernel::copyParametersToContext(ContextImpl& context, const RBTorsionForce& force) {
    if (numTorsions != force.getNumTorsions())
        throw OpenMMException("updateParametersInContext: The number of torsions has changed");

    // Record the values.

    for (int i = 0; i < numTorsions; ++i) {
        int particle1, particle2, particle3, particle4;
        double c0, c1, c2, c3, c4, c5;
        force.getTorsionParameters(i, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
        if (particle1 != torsionIndexArray[i][0] || particle2 != torsionIndexArray[i][1] || particle3 != torsionIndexArray[i][2] || particle4 != torsionIndexArray[i][3])
            throw OpenMMException("updateParametersInContext: The set of particles in a torsion has changed");
peastman's avatar
peastman committed
462
463
464
465
466
467
        torsionParamArray[i][0] = c0;
        torsionParamArray[i][1] = c1;
        torsionParamArray[i][2] = c2;
        torsionParamArray[i][3] = c3;
        torsionParamArray[i][4] = c4;
        torsionParamArray[i][5] = c5;
468
469
470
    }
}

471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
class CpuCalcNonbondedForceKernel::PmeIO : public CalcPmeReciprocalForceKernel::IO {
public:
    PmeIO(float* posq, float* force, int numParticles) : posq(posq), force(force), numParticles(numParticles) {
    }
    float* getPosq() {
        return posq;
    }
    void setForce(float* f) {
        for (int i = 0; i < numParticles; i++) {
            force[4*i] += f[4*i];
            force[4*i+1] += f[4*i+1];
            force[4*i+2] += f[4*i+2];
        }
    }
private:
    float* posq;
    float* force;
    int numParticles;
};

491
CpuNonbondedForce* createCpuNonbondedForceVec(const CpuNeighborList& neighbors);
492
493

CpuCalcNonbondedForceKernel::CpuCalcNonbondedForceKernel(string name, const Platform& platform, CpuPlatform::PlatformData& data) : CalcNonbondedForceKernel(name, platform),
494
        data(data), hasInitializedPme(false), hasInitializedDispersionPme(false), nonbonded(NULL) {
495
496
}

497
CpuCalcNonbondedForceKernel::~CpuCalcNonbondedForceKernel() {
498
499
    if (nonbonded != NULL)
        delete nonbonded;
500
501
502
}

void CpuCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
503
504
    chargePosqIndex = data.requestPosqIndex();
    ljPosqIndex = data.requestPosqIndex();
505
506
507

    // Identify which exceptions are 1-4 interactions.

508
509
510
511
512
513
514
515
    set<int> exceptionsWithOffsets;
    for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
        string param;
        int exception;
        double charge, sigma, epsilon;
        force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
        exceptionsWithOffsets.insert(exception);
    }
516
517
518
519
520
521
522
523
524
    numParticles = force.getNumParticles();
    exclusions.resize(numParticles);
    vector<int> nb14s;
    for (int i = 0; i < force.getNumExceptions(); i++) {
        int particle1, particle2;
        double chargeProd, sigma, epsilon;
        force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
        exclusions[particle1].insert(particle2);
        exclusions[particle2].insert(particle1);
525
526
        if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end()) {
            nb14Index[i] = nb14s.size();
527
            nb14s.push_back(i);
528
        }
529
530
531
532
533
    }

    // Record the particle parameters.

    num14 = nb14s.size();
534
535
    bonded14IndexArray.resize(num14, vector<int>(2));
    bonded14ParamArray.resize(num14, vector<double>(3));
peastman's avatar
peastman committed
536
    particleParams.resize(numParticles);
537
    charges.resize(numParticles);
538
    C6params.resize(numParticles);
539
540
541
    baseParticleParams.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
       force.getParticleParameters(i, baseParticleParams[i][0], baseParticleParams[i][1], baseParticleParams[i][2]);
542
    
543
    // Record exception parameters.
544
    
545
    baseExceptionParams.resize(num14);
546
547
    for (int i = 0; i < num14; ++i) {
        int particle1, particle2;
548
        force.getExceptionParameters(nb14s[i], particle1, particle2, baseExceptionParams[i][0], baseExceptionParams[i][1], baseExceptionParams[i][2]);
549
550
        bonded14IndexArray[i][0] = particle1;
        bonded14IndexArray[i][1] = particle2;
551
    }
peastman's avatar
peastman committed
552
    bondForce.initialize(system.getNumParticles(), num14, 2, bonded14IndexArray, data.threads);
553
    
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
    // Record information about parameter offsets.
    
    hasParticleOffsets = (force.getNumParticleParameterOffsets() > 0);
    hasExceptionOffsets = (force.getNumExceptionParameterOffsets() > 0);
    particleParamOffsets.resize(force.getNumParticles());
    exceptionParamOffsets.resize(force.getNumExceptions());
    for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
        string param;
        int particle;
        double charge, sigma, epsilon;
        force.getParticleParameterOffset(i, param, particle, charge, sigma, epsilon);
        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
        int paramIndex;
        if (paramPos == paramNames.end()) {
            paramIndex = paramNames.size();
            paramNames.push_back(param);
        }
        else
            paramIndex = paramPos-paramNames.begin();
        particleParamOffsets[particle].push_back(make_tuple(charge, sigma, epsilon, paramIndex));
    }
    for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
        string param;
        int exception;
        double charge, sigma, epsilon;
        force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
        int paramIndex;
        if (paramPos == paramNames.end()) {
            paramIndex = paramNames.size();
            paramNames.push_back(param);
        }
        else
            paramIndex = paramPos-paramNames.begin();
588
        exceptionParamOffsets[nb14Index[exception]].push_back(make_tuple(charge, sigma, epsilon, paramIndex));
589
    }
590
    paramValues.resize(paramNames.size(), 0.0);
591

592
593
    // Record other parameters.
    
594
595
    nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
    nonbondedCutoff = force.getCutoffDistance();
596
597
    if (nonbondedMethod == NoCutoff) {
        data.requestNeighborList(0.0, 0.0, true, exclusions);
598
        useSwitchingFunction = false;
599
    }
600
    else {
601
        data.requestNeighborList(nonbondedCutoff, 0.25*nonbondedCutoff, true, exclusions);
602
603
604
605
606
607
608
609
610
611
        useSwitchingFunction = force.getUseSwitchingFunction();
        switchingDistance = force.getSwitchingDistance();
    }
    if (nonbondedMethod == Ewald) {
        double alpha;
        NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmax[0], kmax[1], kmax[2]);
        ewaldAlpha = alpha;
    }
    else if (nonbondedMethod == PME) {
        double alpha;
612
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2], false);
613
614
        ewaldAlpha = alpha;
    }
615
616
617
    else if (nonbondedMethod == LJPME) {
        double alpha;
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2], false);
peastman's avatar
peastman committed
618
        ewaldAlpha = alpha;
619
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, dispersionGridSize[0], dispersionGridSize[1], dispersionGridSize[2], true);
peastman's avatar
peastman committed
620
        ewaldDispersionAlpha = alpha;
621
622
        useSwitchingFunction = false;
    }
623
624
625
626
    if (nonbondedMethod == NoCutoff || nonbondedMethod == CutoffNonPeriodic)
        exceptionsArePeriodic = false;
    else
        exceptionsArePeriodic = force.getExceptionsUsePeriodicBoundaryConditions();
627
628
629
630
631
    rfDielectric = force.getReactionFieldDielectric();
    if (force.getUseDispersionCorrection())
        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
    else
        dispersionCoefficient = 0.0;
632
    data.isPeriodic |= (nonbondedMethod == CutoffPeriodic || nonbondedMethod == Ewald || nonbondedMethod == PME || nonbondedMethod == LJPME);
633
    nonbonded = createCpuNonbondedForceVec(*data.neighborList);
634
635
636
}

double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
637
638
639
    if (!hasInitializedPme) {
        hasInitializedPme = true;
        useOptimizedPme = false;
640
        computeParameters(context, false);
641
642
643
        if (nonbondedMethod == PME) {
            // If available, use the optimized PME implementation.

peastman's avatar
peastman committed
644
645
646
647
            vector<string> kernelNames;
            kernelNames.push_back("CalcPmeReciprocalForce");
            useOptimizedPme = getPlatform().supportsKernels(kernelNames);
            if (useOptimizedPme) {
648
                optimizedPme = getPlatform().createKernel(CalcPmeReciprocalForceKernel::Name(), context);
649
                optimizedPme.getAs<CalcPmeReciprocalForceKernel>().initialize(gridSize[0], gridSize[1], gridSize[2], numParticles, ewaldAlpha, data.deterministicForces);
650
651
            }
        }
652
653
654
655
656
        if (nonbondedMethod == LJPME) {
            // If available, use the optimized PME implementation.

            vector<string> kernelNames;
            kernelNames.push_back("CalcPmeReciprocalForce");
657
            kernelNames.push_back("CalcDispersionPmeReciprocalForce");
658
659
660
            useOptimizedPme = getPlatform().supportsKernels(kernelNames);
            if (useOptimizedPme) {
                optimizedPme = getPlatform().createKernel(CalcPmeReciprocalForceKernel::Name(), context);
661
                optimizedPme.getAs<CalcPmeReciprocalForceKernel>().initialize(gridSize[0], gridSize[1], gridSize[2], numParticles, ewaldAlpha, data.deterministicForces);
662
663
                optimizedDispersionPme = getPlatform().createKernel(CalcDispersionPmeReciprocalForceKernel::Name(), context);
                optimizedDispersionPme.getAs<CalcDispersionPmeReciprocalForceKernel>().initialize(dispersionGridSize[0], dispersionGridSize[1],
664
                                                                                                  dispersionGridSize[2], numParticles, ewaldDispersionAlpha, data.deterministicForces);
665
666
            }
        }
667
    }
668
    computeParameters(context, true);
669
    copyChargesToPosq(context, charges, chargePosqIndex);
670
    AlignedArray<float>& posq = data.posq;
peastman's avatar
peastman committed
671
672
673
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    Vec3* boxVectors = extractBoxVectors(context);
674
    double energy = (includeReciprocal ? ewaldSelfEnergy : 0.0);
675
676
    bool ewald  = (nonbondedMethod == Ewald);
    bool pme  = (nonbondedMethod == PME);
677
    bool ljpme = (nonbondedMethod == LJPME);
678
    if (nonbondedMethod != NoCutoff)
679
        nonbonded->setUseCutoff(nonbondedCutoff, rfDielectric);
680
    if (data.isPeriodic) {
peastman's avatar
peastman committed
681
        Vec3* boxVectors = extractBoxVectors(context);
682
        double minAllowedSize = 1.999999*nonbondedCutoff;
683
        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
684
            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
685
        nonbonded->setPeriodic(boxVectors);
686
        nonbonded->setPeriodicExceptions(exceptionsArePeriodic);
687
688
    }
    if (ewald)
689
        nonbonded->setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
690
    if (pme)
691
        nonbonded->setUsePME(ewaldAlpha, gridSize);
692
    if (useSwitchingFunction)
693
        nonbonded->setUseSwitchingFunction(switchingDistance);
694
695
696
697
    if (ljpme){
        nonbonded->setUsePME(ewaldAlpha, gridSize);
        nonbonded->setUseLJPME(ewaldDispersionAlpha, dispersionGridSize);
    }
698
    double nonbondedEnergy = 0;
peastman's avatar
peastman committed
699
    if (includeDirect)
700
        nonbonded->calculateDirectIxn(numParticles, &posq[0], posData, particleParams, C6params, exclusions, data.threadForce, includeEnergy ? &nonbondedEnergy : NULL, data.threads);
701
702
    if (includeReciprocal) {
        if (useOptimizedPme) {
703
            PmeIO io(&posq[0], &data.threadForce[0][0], numParticles);
peastman's avatar
peastman committed
704
705
            Vec3 periodicBoxVectors[3] = {boxVectors[0], boxVectors[1], boxVectors[2]};
            optimizedPme.getAs<CalcPmeReciprocalForceKernel>().beginComputation(io, periodicBoxVectors, includeEnergy);
706
            nonbondedEnergy += optimizedPme.getAs<CalcPmeReciprocalForceKernel>().finishComputation(io);
707
708
709
710
711
            if (nonbondedMethod == LJPME) {
                copyChargesToPosq(context, C6params, ljPosqIndex);
                optimizedDispersionPme.getAs<CalcDispersionPmeReciprocalForceKernel>().beginComputation(io, periodicBoxVectors, includeEnergy);
                nonbondedEnergy += optimizedDispersionPme.getAs<CalcDispersionPmeReciprocalForceKernel>().finishComputation(io);
            }
712
713
        }
        else
714
            nonbonded->calculateReciprocalIxn(numParticles, &posq[0], posData, particleParams, C6params, exclusions, forceData, includeEnergy ? &nonbondedEnergy : NULL);
715
    }
peastman's avatar
peastman committed
716
    energy += nonbondedEnergy;
717
718
    if (includeDirect) {
        ReferenceLJCoulomb14 nonbonded14;
719
720
721
722
        if (exceptionsArePeriodic) {
            Vec3* boxVectors = extractBoxVectors(context);
            nonbonded14.setPeriodic(boxVectors);
        }
peastman's avatar
peastman committed
723
        bondForce.calculateForce(posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
724
        if (data.isPeriodic && nonbondedMethod != LJPME)
725
            energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
726
727
    }
    return energy;
728
729
}

730
void CpuCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const NonbondedForce& force, int firstParticle, int lastParticle, int firstException, int lastException) {
731
732
    if (force.getNumParticles() != numParticles)
        throw OpenMMException("updateParametersInContext: The number of particles has changed");
peastman's avatar
peastman committed
733
734
735
736
737
738
739
740
741
742
743

    // Identify which exceptions are 1-4 interactions.

    set<int> exceptionsWithOffsets;
    for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
        string param;
        int exception;
        double charge, sigma, epsilon;
        force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
        exceptionsWithOffsets.insert(exception);
    }
744
745
746
747
748
    vector<int> nb14s;
    for (int i = 0; i < force.getNumExceptions(); i++) {
        int particle1, particle2;
        double chargeProd, sigma, epsilon;
        force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
749
750
751
752
753
        if (nb14Index.find(i) == nb14Index.end()) {
            if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end())
                throw OpenMMException("updateParametersInContext: The set of non-excluded exceptions has changed");
        }
        else
754
755
756
757
758
759
760
            nb14s.push_back(i);
    }
    if (nb14s.size() != num14)
        throw OpenMMException("updateParametersInContext: The number of non-excluded exceptions has changed");

    // Record the values.

761
    for (int i = firstParticle; i <= lastParticle; ++i)
762
       force.getParticleParameters(i, baseParticleParams[i][0], baseParticleParams[i][1], baseParticleParams[i][2]);
763
764
    for (int i = 0; i < num14; ++i) {
        int particle1, particle2;
765
        force.getExceptionParameters(nb14s[i], particle1, particle2, baseExceptionParams[i][0], baseExceptionParams[i][1], baseExceptionParams[i][2]);
766
767
768
        bonded14IndexArray[i][0] = particle1;
        bonded14IndexArray[i][1] = particle2;
    }
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
    particleParamOffsets.clear();
    exceptionParamOffsets.clear();
    particleParamOffsets.resize(force.getNumParticles());
    exceptionParamOffsets.resize(force.getNumExceptions());
    for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
        string param;
        int particle;
        double charge, sigma, epsilon;
        force.getParticleParameterOffset(i, param, particle, charge, sigma, epsilon);
        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
        if (paramPos == paramNames.end())
            throw OpenMMException("updateParametersInContext: The parameter of a particle parameter offset has changed");
        int paramIndex = paramPos-paramNames.begin();
        particleParamOffsets[particle].push_back(make_tuple(charge, sigma, epsilon, paramIndex));
    }
    for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
        string param;
        int exception;
        double charge, sigma, epsilon;
        force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
        auto paramPos = find(paramNames.begin(), paramNames.end(), param);
        if (paramPos == paramNames.end())
            throw OpenMMException("updateParametersInContext: The parameter of an exception parameter offset has changed");
        int paramIndex = paramPos-paramNames.begin();
        exceptionParamOffsets[nb14Index[exception]].push_back(make_tuple(charge, sigma, epsilon, paramIndex));
    }
795
    computeParameters(context, false);
796
797
798
799
800
801
    
    // Recompute the coefficient for the dispersion correction.

    NonbondedForce::NonbondedMethod method = force.getNonbondedMethod();
    if (force.getUseDispersionCorrection() && (method == NonbondedForce::CutoffPeriodic || method == NonbondedForce::Ewald || method == NonbondedForce::PME))
        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
802
}
803

804
void CpuCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
805
    if (nonbondedMethod != PME && nonbondedMethod != LJPME)
806
807
808
809
810
811
812
813
814
815
816
        throw OpenMMException("getPMEParametersInContext: This Context is not using PME");
    if (useOptimizedPme)
        optimizedPme.getAs<const CalcPmeReciprocalForceKernel>().getPMEParameters(alpha, nx, ny, nz);
    else {
        alpha = ewaldAlpha;
        nx = gridSize[0];
        ny = gridSize[1];
        nz = gridSize[2];
    }
}

817
void CpuCalcNonbondedForceKernel::getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
818
819
820
    if (nonbondedMethod != LJPME)
        throw OpenMMException("getPMEParametersInContext: This Context is not using PME");
    if (useOptimizedPme)
821
        optimizedDispersionPme.getAs<const CalcPmeReciprocalForceKernel>().getPMEParameters(alpha, nx, ny, nz);
822
    else {
823
824
825
826
        alpha = ewaldDispersionAlpha;
        nx = dispersionGridSize[0];
        ny = dispersionGridSize[1];
        nz = dispersionGridSize[2];
827
828
829
    }
}

830
void CpuCalcNonbondedForceKernel::computeParameters(ContextImpl& context, bool offsetsOnly) {
831
832
833
834
835
836
837
838
839
840
    bool paramChanged = false;
    for (int i = 0; i < paramNames.size(); i++) {
        double value = context.getParameter(paramNames[i]);
        if (value != paramValues[i]) {
            paramValues[i] = value;;
            paramChanged = true;
        }
    }
    if (!paramChanged && offsetsOnly)
        return;
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868

    // Compute particle parameters.

    if (hasParticleOffsets || !offsetsOnly) {
        double sumSquaredCharges = 0.0;
        for (int i = 0; i < numParticles; i++) {
            double charge = baseParticleParams[i][0];
            double sigma = baseParticleParams[i][1];
            double epsilon = baseParticleParams[i][2];
            for (auto& offset : particleParamOffsets[i]) {
                double value = paramValues[get<3>(offset)];
                charge += value*get<0>(offset);
                sigma += value*get<1>(offset);
                epsilon += value*get<2>(offset);
            }
            charges[i] = (float) charge;
            particleParams[i] = make_pair((float) (0.5*sigma), (float) (2.0*sqrt(epsilon)));
            C6params[i] = 8.0*pow(particleParams[i].first, 3.0) * particleParams[i].second;
            sumSquaredCharges += charge*charge;
        }
        if (nonbondedMethod == Ewald || nonbondedMethod == PME || nonbondedMethod == LJPME) {
            ewaldSelfEnergy = -ONE_4PI_EPS0*ewaldAlpha*sumSquaredCharges/sqrt(M_PI);
            if (nonbondedMethod == LJPME)
                for (int atom = 0; atom < numParticles; atom++)
                    ewaldSelfEnergy += pow(ewaldDispersionAlpha, 6.0) * C6params[atom]*C6params[atom] / 12.0;
        }
        else
            ewaldSelfEnergy = 0.0;
869
870
        chargePosqIndex = data.requestPosqIndex();
        ljPosqIndex = data.requestPosqIndex();
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
    }

    // Compute exception parameters.

    if (hasExceptionOffsets || !offsetsOnly) {
        for (int i = 0; i < num14; i++) {
            double chargeProd = baseExceptionParams[i][0];
            double sigma = baseExceptionParams[i][1];
            double epsilon = baseExceptionParams[i][2];
            for (auto& offset : exceptionParamOffsets[i]) {
                double value = paramValues[get<3>(offset)];
                chargeProd += value*get<0>(offset);
                sigma += value*get<1>(offset);
                epsilon += value*get<2>(offset);
            }
            bonded14ParamArray[i][0] = sigma;
            bonded14ParamArray[i][1] = 4.0*epsilon;
            bonded14ParamArray[i][2] = chargeProd;
        }
    }
}

893
CpuCalcCustomNonbondedForceKernel::CpuCalcCustomNonbondedForceKernel(string name, const Platform& platform, CpuPlatform::PlatformData& data) :
894
            CalcCustomNonbondedForceKernel(name, platform), data(data), forceCopy(NULL), nonbonded(NULL) {
895
896
}

897
CpuCalcCustomNonbondedForceKernel::~CpuCalcCustomNonbondedForceKernel() {
898
899
    if (nonbonded != NULL)
        delete nonbonded;
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
    if (forceCopy != NULL)
        delete forceCopy;
}

void CpuCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {

    // Record the exclusions.

    numParticles = force.getNumParticles();
    exclusions.resize(numParticles);
    for (int i = 0; i < force.getNumExclusions(); i++) {
        int particle1, particle2;
        force.getExclusionParticles(i, particle1, particle2);
        exclusions[particle1].insert(particle2);
        exclusions[particle2].insert(particle1);
    }

    // Build the arrays.

919
920
921
    particleParamArray.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        force.getParticleParameters(i, particleParamArray[i]);
922
923
    nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
    nonbondedCutoff = force.getCutoffDistance();
924
925
    if (nonbondedMethod == NoCutoff) {
        data.requestNeighborList(0.0, 0.0, true, exclusions);
926
        useSwitchingFunction = false;
927
    }
928
    else {
929
        data.requestNeighborList(nonbondedCutoff, 0.25*nonbondedCutoff, true, exclusions);
930
931
932
933
        useSwitchingFunction = force.getUseSwitchingFunction();
        switchingDistance = force.getSwitchingDistance();
    }

934
    // Record the tabulated function update counts for future reference.
935
936

    for (int i = 0; i < force.getNumTabulatedFunctions(); i++)
937
        tabulatedFunctionUpdateCount[force.getTabulatedFunctionName(i)] = force.getTabulatedFunction(i).getUpdateCount();
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964

    // Record information for the long range correction.

    if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic && force.getUseLongRangeCorrection()) {
        forceCopy = new CustomNonbondedForce(force);
        hasInitializedLongRangeCorrection = false;
    }
    else {
        longRangeCoefficient = 0.0;
        hasInitializedLongRangeCorrection = true;
    }

    // Record the interaction groups.

    for (int i = 0; i < force.getNumInteractionGroups(); i++) {
        set<int> set1, set2;
        force.getInteractionGroupParameters(i, set1, set2);
        interactionGroups.push_back(make_pair(set1, set2));
    }
    data.isPeriodic |= (nonbondedMethod == CutoffPeriodic);

    // Create the interaction.

    createInteraction(force);
}

void CpuCalcCustomNonbondedForceKernel::createInteraction(const CustomNonbondedForce& force) {
965
966
967
    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
968
    for (int i = 0; i < force.getNumTabulatedFunctions(); i++)
969
970
971
972
        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));

    // Parse the various expressions used to calculate the force.

973
974
    Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
    Lepton::ParsedExpression forceExpression = energyExpression.differentiate("r");
975
    for (int i = 0; i < force.getNumPerParticleParameters(); i++)
976
977
978
979
980
        parameterNames.push_back(force.getPerParticleParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
        globalParameterNames.push_back(force.getGlobalParameterName(i));
        globalParamValues[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
    }
981
982
983
984
985
986
987
988
989
990
    set<string> particleVariables, pairVariables;
    particleVariables.insert("r");
    pairVariables.insert("r");
    for (auto& name : parameterNames) {
        particleVariables.insert(name);
        pairVariables.insert(name+"1");
        pairVariables.insert(name+"2");
    }
    particleVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
    pairVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
991
    vector<Lepton::ParsedExpression> computedValueExpressions, energyParamDerivExpressions;
992
993
994
995
996
997
    for (int i = 0; i < force.getNumComputedValues(); i++) {
        string name, exp;
        force.getComputedValueParameters(i, name, exp);
        Lepton::ParsedExpression parsed = Lepton::Parser::parse(exp, functions);
        validateVariables(parsed.getRootNode(), particleVariables);
        computedValueNames.push_back(name);
998
        computedValueExpressions.push_back(parsed);
999
    }
1000
1001
1002
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
        string param = force.getEnergyParameterDerivativeName(i);
        energyParamDerivNames.push_back(param);
1003
        energyParamDerivExpressions.push_back(energyExpression.differentiate(param));
1004
    }
1005
1006
1007
    for (auto& name : computedValueNames) {
        pairVariables.insert(name+"1");
        pairVariables.insert(name+"2");
1008
    }
1009
    validateVariables(energyExpression.getRootNode(), pairVariables);
1010
1011
1012

    // Delete the custom functions.

peastman's avatar
peastman committed
1013
1014
    for (auto& function : functions)
        delete function.second;
1015
1016
1017

    // Create the object that computes the interaction.

1018
    nonbonded = createCpuCustomNonbondedForce(data.threads, *data.neighborList);
1019
1020
    nonbonded->initialize(energyExpression, forceExpression, parameterNames, exclusions, energyParamDerivExpressions,
            computedValueNames, computedValueExpressions);
peastman's avatar
Bug fix  
peastman committed
1021
1022
    if (interactionGroups.size() > 0)
        nonbonded->setInteractionGroups(interactionGroups);
1023
1024
1025
}

double CpuCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1026
1027
1028
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    Vec3* boxVectors = extractBoxVectors(context);
1029
1030
    double energy = 0;
    bool periodic = (nonbondedMethod == CutoffPeriodic);
1031
    if (nonbondedMethod != NoCutoff)
1032
        nonbonded->setUseCutoff(nonbondedCutoff);
1033
1034
    if (periodic) {
        double minAllowedSize = 2*nonbondedCutoff;
1035
        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
1036
            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
1037
        nonbonded->setPeriodic(boxVectors);
1038
1039
    }
    bool globalParamsChanged = false;
peastman's avatar
peastman committed
1040
1041
1042
    for (auto& name : globalParameterNames) {
        double value = context.getParameter(name);
        if (globalParamValues[name] != value)
1043
            globalParamsChanged = true;
peastman's avatar
peastman committed
1044
        globalParamValues[name] = value;
1045
1046
    }
    if (useSwitchingFunction)
1047
        nonbonded->setUseSwitchingFunction(switchingDistance);
1048
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
1049
    nonbonded->calculatePairIxn(numParticles, &data.posq[0], posData, particleParamArray, globalParamValues, data.threadForce, includeForces, includeEnergy, energy, &energyParamDerivValues[0]);
1050
1051
1052
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
1053
1054
1055
    
    // Add in the long range correction.
    
1056
    if (!hasInitializedLongRangeCorrection) {
1057
        longRangeCorrectionData = CustomNonbondedForceImpl::prepareLongRangeCorrection(*forceCopy, data.threads.getNumThreads());
1058
        CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, longRangeCorrectionData, context.getOwner(), longRangeCoefficient, longRangeCoefficientDerivs, data.threads);
1059
1060
        hasInitializedLongRangeCorrection = true;
    }
1061
1062
    else if (globalParamsChanged && forceCopy != NULL)
        CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, longRangeCorrectionData, context.getOwner(), longRangeCoefficient, longRangeCoefficientDerivs, data.threads);
1063
1064
1065
1066
    double volume = boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2];
    energy += longRangeCoefficient/volume;
    for (int i = 0; i < longRangeCoefficientDerivs.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += longRangeCoefficientDerivs[i]/volume;
1067
1068
1069
    return energy;
}

1070
void CpuCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force, int firstParticle, int lastParticle) {
1071
1072
1073
1074
1075
1076
    if (numParticles != force.getNumParticles())
        throw OpenMMException("updateParametersInContext: The number of particles has changed");

    // Record the values.

    int numParameters = force.getNumPerParticleParameters();
1077
1078
    vector<double> parameters;
    for (int i = firstParticle; i <= lastParticle; ++i) {
1079
1080
1081
1082
1083
1084
1085
1086
        force.getParticleParameters(i, parameters);
        for (int j = 0; j < numParameters; j++)
            particleParamArray[i][j] = parameters[j];
    }
    
    // If necessary, recompute the long range correction.
    
    if (forceCopy != NULL) {
1087
        longRangeCorrectionData = CustomNonbondedForceImpl::prepareLongRangeCorrection(force, data.threads.getNumThreads());
1088
        CustomNonbondedForceImpl::calcLongRangeCorrection(force, longRangeCorrectionData, context.getOwner(), longRangeCoefficient, longRangeCoefficientDerivs, data.threads);
1089
1090
1091
        hasInitializedLongRangeCorrection = true;
        *forceCopy = force;
    }
1092
1093
1094
1095
1096
1097

    // See if any tabulated functions have changed.

    bool changed = false;
    for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
        string name = force.getTabulatedFunctionName(i);
1098
1099
        if (force.getTabulatedFunction(i).getUpdateCount() != tabulatedFunctionUpdateCount[name]) {
            tabulatedFunctionUpdateCount[name] = force.getTabulatedFunction(i).getUpdateCount();
1100
1101
1102
1103
1104
1105
1106
1107
            changed = true;
        }
    }
    if (changed) {
        delete nonbonded;
        nonbonded = NULL;
        createInteraction(force);
    }
1108
1109
}

1110
1111
1112
1113
CpuCalcGBSAOBCForceKernel::~CpuCalcGBSAOBCForceKernel() {
}

void CpuCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
1114
    posqIndex = data.requestPosqIndex();
1115
1116
    int numParticles = system.getNumParticles();
    particleParams.resize(numParticles);
1117
    charges.resize(numParticles);
1118
1119
1120
    for (int i = 0; i < numParticles; ++i) {
        double charge, radius, scalingFactor;
        force.getParticleParameters(i, charge, radius, scalingFactor);
1121
        charges[i] = (float) charge;
1122
1123
        radius -= 0.009;
        particleParams[i] = make_pair((float) radius, (float) (scalingFactor*radius));
1124
1125
1126
1127
    }
    obc.setParticleParameters(particleParams);
    obc.setSolventDielectric((float) force.getSolventDielectric());
    obc.setSoluteDielectric((float) force.getSoluteDielectric());
1128
    obc.setSurfaceAreaEnergy((float) force.getSurfaceAreaEnergy());
1129
1130
    if (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff)
        obc.setUseCutoff((float) force.getCutoffDistance());
1131
    data.isPeriodic |= (force.getNonbondedMethod() == GBSAOBCForce::CutoffPeriodic);
1132
1133
1134
}

double CpuCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
1135
    copyChargesToPosq(context, charges, posqIndex);
1136
    if (data.isPeriodic) {
peastman's avatar
peastman committed
1137
        Vec3& boxSize = extractBoxSize(context);
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
        float floatBoxSize[3] = {(float) boxSize[0], (float) boxSize[1], (float) boxSize[2]};
        obc.setPeriodic(floatBoxSize);
    }
    double energy = 0.0;
    obc.computeForce(data.posq, data.threadForce, includeEnergy ? &energy : NULL, data.threads);
    return energy;
}

void CpuCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force) {
    int numParticles = force.getNumParticles();
    if (numParticles != obc.getParticleParameters().size())
        throw OpenMMException("updateParametersInContext: The number of particles has changed");

    // Record the values.

1153
    posqIndex = data.requestPosqIndex();
1154
1155
1156
    for (int i = 0; i < numParticles; ++i) {
        double charge, radius, scalingFactor;
        force.getParticleParameters(i, charge, radius, scalingFactor);
1157
        charges[i] = (float) charge;
1158
1159
        radius -= 0.009;
        particleParams[i] = make_pair((float) radius, (float) (scalingFactor*radius));
1160
1161
1162
    }
    obc.setParticleParameters(particleParams);
}
1163

1164
CpuCalcCustomGBForceKernel::~CpuCalcCustomGBForceKernel() {
1165
1166
    if (ixn != NULL)
        delete ixn;
Peter Eastman's avatar
Peter Eastman committed
1167
1168
    if (neighborList != NULL)
        delete neighborList;
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
}

void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
    if (force.getNumComputedValues() > 0) {
        string name, expression;
        CustomGBForce::ComputationType type;
        force.getComputedValueParameters(0, name, expression, type);
        if (type == CustomGBForce::SingleParticle)
            throw OpenMMException("CpuPlatform requires that the first computed value for a CustomGBForce be of type ParticlePair or ParticlePairNoExclusions.");
        for (int i = 1; i < force.getNumComputedValues(); i++) {
            force.getComputedValueParameters(i, name, expression, type);
            if (type != CustomGBForce::SingleParticle)
                throw OpenMMException("CpuPlatform requires that a CustomGBForce only have one computed value of type ParticlePair or ParticlePairNoExclusions.");
        }
    }

    // Record the exclusions.

    numParticles = force.getNumParticles();
    exclusions.resize(numParticles);
    for (int i = 0; i < force.getNumExclusions(); i++) {
        int particle1, particle2;
        force.getExclusionParticles(i, particle1, particle2);
        exclusions[particle1].insert(particle2);
        exclusions[particle2].insert(particle1);
    }

    // Build the arrays.

1198
1199
1200
    particleParamArray.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        force.getParticleParameters(i, particleParamArray[i]);
1201
    for (int i = 0; i < force.getNumPerParticleParameters(); i++)
1202
1203
1204
1205
        particleParameterNames.push_back(force.getPerParticleParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
    nonbondedMethod = CalcCustomGBForceKernel::NonbondedMethod(force.getNonbondedMethod());
peastman's avatar
peastman committed
1206
    nonbondedCutoff = force.getCutoffDistance();
1207
    if (nonbondedMethod != NoCutoff)
Peter Eastman's avatar
Peter Eastman committed
1208
        neighborList = new CpuNeighborList(4);
1209
1210
    data.isPeriodic |= (force.getNonbondedMethod() == CustomGBForce::CutoffPeriodic);

1211
    // Record the tabulated function update counts for future reference.
1212
1213

    for (int i = 0; i < force.getNumTabulatedFunctions(); i++)
1214
        tabulatedFunctionUpdateCount[force.getTabulatedFunctionName(i)] = force.getTabulatedFunction(i).getUpdateCount();
1215
1216

    // Create the interaction.
1217

1218
1219
1220
1221
    createInteraction(force);
}

void CpuCalcCustomGBForceKernel::createInteraction(const CustomGBForce& force) {
1222
1223
1224
    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
1225
    for (int i = 0; i < force.getNumTabulatedFunctions(); i++)
1226
1227
1228
1229
        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));

    // Parse the expressions for computed values.

1230
1231
1232
    valueTypes.clear();
    valueNames.clear();
    energyParamDerivNames.clear();
1233
1234
    vector<vector<Lepton::CompiledExpression> > valueDerivExpressions(force.getNumComputedValues());
    vector<vector<Lepton::CompiledExpression> > valueGradientExpressions(force.getNumComputedValues());
1235
    vector<vector<Lepton::CompiledExpression> > valueParamDerivExpressions(force.getNumComputedValues());
1236
1237
    vector<Lepton::CompiledExpression> valueExpressions;
    vector<Lepton::CompiledExpression> energyExpressions;
1238
1239
1240
1241
1242
    set<string> particleVariables, pairVariables;
    pairVariables.insert("r");
    particleVariables.insert("x");
    particleVariables.insert("y");
    particleVariables.insert("z");
1243
    for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
1244
1245
1246
1247
1248
1249
        particleVariables.insert(particleParameterNames[i]);
        pairVariables.insert(particleParameterNames[i]+"1");
        pairVariables.insert(particleParameterNames[i]+"2");
    }
    particleVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
    pairVariables.insert(globalParameterNames.begin(), globalParameterNames.end());
1250
1251
1252
1253
1254
1255
1256
1257
    for (int i = 0; i < force.getNumComputedValues(); i++) {
        string name, expression;
        CustomGBForce::ComputationType type;
        force.getComputedValueParameters(i, name, expression, type);
        Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
        valueExpressions.push_back(ex.createCompiledExpression());
        valueTypes.push_back(type);
        valueNames.push_back(name);
1258
        if (i == 0) {
1259
            valueDerivExpressions[i].push_back(ex.differentiate("r").createCompiledExpression());
1260
1261
            validateVariables(ex.getRootNode(), pairVariables);
        }
1262
1263
1264
1265
1266
1267
        else {
            valueGradientExpressions[i].push_back(ex.differentiate("x").createCompiledExpression());
            valueGradientExpressions[i].push_back(ex.differentiate("y").createCompiledExpression());
            valueGradientExpressions[i].push_back(ex.differentiate("z").createCompiledExpression());
            for (int j = 0; j < i; j++)
                valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).createCompiledExpression());
1268
            validateVariables(ex.getRootNode(), particleVariables);
1269
        }
1270
1271
1272
1273
1274
        for (int j = 0; j < force.getNumEnergyParameterDerivatives(); j++) {
            string param = force.getEnergyParameterDerivativeName(j);
            energyParamDerivNames.push_back(param);
            valueParamDerivExpressions[i].push_back(ex.differentiate(param).createCompiledExpression());
        }
1275
1276
1277
        particleVariables.insert(name);
        pairVariables.insert(name+"1");
        pairVariables.insert(name+"2");
1278
1279
1280
1281
    }

    // Parse the expressions for energy terms.

1282
    energyTypes.clear();
1283
1284
    vector<vector<Lepton::CompiledExpression> > energyDerivExpressions(force.getNumEnergyTerms());
    vector<vector<Lepton::CompiledExpression> > energyGradientExpressions(force.getNumEnergyTerms());
1285
    vector<vector<Lepton::CompiledExpression> > energyParamDerivExpressions(force.getNumEnergyTerms());
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
    for (int i = 0; i < force.getNumEnergyTerms(); i++) {
        string expression;
        CustomGBForce::ComputationType type;
        force.getEnergyTermParameters(i, expression, type);
        Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
        energyExpressions.push_back(ex.createCompiledExpression());
        energyTypes.push_back(type);
        if (type != CustomGBForce::SingleParticle)
            energyDerivExpressions[i].push_back(ex.differentiate("r").createCompiledExpression());
        for (int j = 0; j < force.getNumComputedValues(); j++) {
            if (type == CustomGBForce::SingleParticle) {
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).createCompiledExpression());
                energyGradientExpressions[i].push_back(ex.differentiate("x").createCompiledExpression());
                energyGradientExpressions[i].push_back(ex.differentiate("y").createCompiledExpression());
                energyGradientExpressions[i].push_back(ex.differentiate("z").createCompiledExpression());
1301
                validateVariables(ex.getRootNode(), particleVariables);
1302
1303
1304
1305
            }
            else {
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").createCompiledExpression());
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").createCompiledExpression());
1306
                validateVariables(ex.getRootNode(), pairVariables);
1307
1308
            }
        }
1309
1310
        for (int j = 0; j < force.getNumEnergyParameterDerivatives(); j++)
            energyParamDerivExpressions[i].push_back(ex.differentiate(force.getEnergyParameterDerivativeName(j)).createCompiledExpression());
1311
1312
1313
1314
    }

    // Delete the custom functions.

peastman's avatar
peastman committed
1315
1316
    for (auto& function : functions)
        delete function.second;
1317
1318
1319
    ixn = new CpuCustomGBForce(numParticles, exclusions, valueExpressions, valueDerivExpressions, valueGradientExpressions, valueParamDerivExpressions,
        valueNames, valueTypes, energyExpressions, energyDerivExpressions, energyGradientExpressions, energyParamDerivExpressions, energyTypes,
        particleParameterNames, data.threads);
1320
1321
1322
}

double CpuCalcCustomGBForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1323
1324
1325
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
    Vec3* boxVectors = extractBoxVectors(context);
1326
1327
    if (data.isPeriodic)
        ixn->setPeriodic(extractBoxSize(context));
1328
    if (nonbondedMethod != NoCutoff) {
1329
        vector<set<int> > noExclusions(numParticles);
Peter Eastman's avatar
Peter Eastman committed
1330
1331
        neighborList->computeNeighborList(numParticles, data.posq, noExclusions, boxVectors, data.isPeriodic, nonbondedCutoff, data.threads);
        ixn->setUseCutoff(nonbondedCutoff, *neighborList);
1332
1333
    }
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1334
1335
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1336
1337
1338
1339
1340
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
    ixn->calculateIxn(numParticles, &data.posq[0], particleParamArray, globalParameters, data.threadForce, includeForces, includeEnergy, energy, &energyParamDerivValues[0]);
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
    return energy;
}

void CpuCalcCustomGBForceKernel::copyParametersToContext(ContextImpl& context, const CustomGBForce& force) {
    if (numParticles != force.getNumParticles())
        throw OpenMMException("updateParametersInContext: The number of particles has changed");

    // Record the values.

    int numParameters = force.getNumPerParticleParameters();
    vector<double> params;
    for (int i = 0; i < numParticles; ++i) {
        vector<double> parameters;
        force.getParticleParameters(i, parameters);
        for (int j = 0; j < numParameters; j++)
peastman's avatar
peastman committed
1356
            particleParamArray[i][j] = static_cast<double>(parameters[j]);
1357
    }
1358
1359
1360
1361
1362
1363

    // See if any tabulated functions have changed.

    bool changed = false;
    for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
        string name = force.getTabulatedFunctionName(i);
1364
1365
        if (force.getTabulatedFunction(i).getUpdateCount() != tabulatedFunctionUpdateCount[name]) {
            tabulatedFunctionUpdateCount[name] = force.getTabulatedFunction(i).getUpdateCount();
1366
1367
1368
1369
1370
1371
1372
1373
            changed = true;
        }
    }
    if (changed) {
        delete ixn;
        ixn = NULL;
        createInteraction(force);
    }
1374
1375
}

1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
CpuCalcCustomManyParticleForceKernel::~CpuCalcCustomManyParticleForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

void CpuCalcCustomManyParticleForceKernel::initialize(const System& system, const CustomManyParticleForce& force) {

    // Build the arrays.

    numParticles = system.getNumParticles();
1386
    particleParamArray.resize(numParticles);
1387
1388
    for (int i = 0; i < numParticles; ++i) {
        int type;
1389
        force.getParticleParameters(i, particleParamArray[i], type);
1390
1391
1392
    }
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
1393

1394
    // Record the tabulated function update counts for future reference.
1395
1396

    for (int i = 0; i < force.getNumTabulatedFunctions(); i++)
1397
        tabulatedFunctionUpdateCount[force.getTabulatedFunctionName(i)] = force.getTabulatedFunction(i).getUpdateCount();
1398
1399
1400

    // Create the interaction.

1401
    ixn = new CpuCustomManyParticleForce(force, data.threads);
1402
1403
    nonbondedMethod = CalcCustomManyParticleForceKernel::NonbondedMethod(force.getNonbondedMethod());
    cutoffDistance = force.getCutoffDistance();
1404
    data.isPeriodic |= (nonbondedMethod == CutoffPeriodic);
1405
1406
1407
1408
}

double CpuCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1409
1410
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1411
    if (nonbondedMethod == CutoffPeriodic) {
peastman's avatar
peastman committed
1412
        Vec3* boxVectors = extractBoxVectors(context);
1413
        double minAllowedSize = 2*cutoffDistance;
1414
        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
1415
            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
1416
        ixn->setPeriodic(boxVectors);
1417
    }
1418
    double energy = 0;
1419
    ixn->calculateIxn(data.posq, particleParamArray, globalParameters, data.threadForce, includeForces, includeEnergy, energy);
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
    return energy;
}

void CpuCalcCustomManyParticleForceKernel::copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force) {
    if (numParticles != force.getNumParticles())
        throw OpenMMException("updateParametersInContext: The number of particles has changed");

    // Record the values.

    int numParameters = force.getNumPerParticleParameters();
    vector<double> params;
    for (int i = 0; i < numParticles; ++i) {
        vector<double> parameters;
        int type;
        force.getParticleParameters(i, parameters, type);
        for (int j = 0; j < numParameters; j++)
peastman's avatar
peastman committed
1436
            particleParamArray[i][j] = static_cast<double>(parameters[j]);
1437
    }
1438
1439
1440
1441
1442
1443

    // See if any tabulated functions have changed.

    bool changed = false;
    for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
        string name = force.getTabulatedFunctionName(i);
1444
1445
        if (force.getTabulatedFunction(i).getUpdateCount() != tabulatedFunctionUpdateCount[name]) {
            tabulatedFunctionUpdateCount[name] = force.getTabulatedFunction(i).getUpdateCount();
1446
1447
1448
1449
1450
1451
1452
1453
            changed = true;
        }
    }
    if (changed) {
        delete ixn;
        ixn = NULL;
        ixn = new CpuCustomManyParticleForce(force, data.threads);
    }
1454
1455
}

1456
1457
1458
1459
1460
1461
1462
CpuCalcGayBerneForceKernel::~CpuCalcGayBerneForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

void CpuCalcGayBerneForceKernel::initialize(const System& system, const GayBerneForce& force) {
    ixn = new CpuGayBerneForce(force);
1463
    data.isPeriodic |= (force.getNonbondedMethod() == GayBerneForce::CutoffPeriodic);
1464
1465
1466
1467
    if (force.getNonbondedMethod() != GayBerneForce::NoCutoff) {
        double cutoff = force.getCutoffDistance();
        data.requestNeighborList(cutoff, 0.1*cutoff, true, ixn->getExclusions());
    }
1468
1469
1470
}

double CpuCalcGayBerneForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
1471
    return ixn->calculateForce(extractPositions(context), extractForces(context), data.threadForce, extractBoxVectors(context), data);
1472
1473
1474
1475
1476
1477
1478
1479
}

void CpuCalcGayBerneForceKernel::copyParametersToContext(ContextImpl& context, const GayBerneForce& force) {
    delete ixn;
    ixn = NULL;
    ixn = new CpuGayBerneForce(force);
}

1480
CpuIntegrateLangevinMiddleStepKernel::~CpuIntegrateLangevinMiddleStepKernel() {
1481
1482
1483
1484
    if (dynamics)
        delete dynamics;
}

1485
void CpuIntegrateLangevinMiddleStepKernel::initialize(const System& system, const LangevinMiddleIntegrator& integrator) {
1486
1487
1488
1489
1490
1491
1492
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        masses[i] = system.getParticleMass(i);
    data.random.initialize(integrator.getRandomNumberSeed(), data.threads.getNumThreads());
}

1493
void CpuIntegrateLangevinMiddleStepKernel::execute(ContextImpl& context, const LangevinMiddleIntegrator& integrator) {
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics)
            delete dynamics;
1504
        dynamics = new CpuLangevinMiddleDynamics(context.getSystem().getNumParticles(), stepSize, friction, temperature, data.threads, data.random);
1505
        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
1506
        dynamics->setVirtualSites(extractVirtualSites(context));
1507
1508
1509
1510
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
1511
    dynamics->update(context, posData, velData, masses, integrator.getConstraintTolerance());
1512
1513
1514
1515
1516
    ReferencePlatform::PlatformData* refData = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
    refData->time += stepSize;
    refData->stepCount++;
}

1517
double CpuIntegrateLangevinMiddleStepKernel::computeKineticEnergy(ContextImpl& context, const LangevinMiddleIntegrator& integrator) {
1518
1519
    return computeShiftedKineticEnergy(context, masses, 0.0);
}