ReferenceKernels.cpp 108 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) 2008-2018 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 "ReferenceKernels.h"
33
#include "ReferenceObc.h"
34
35
36
37
38
39
#include "ReferenceAndersenThermostat.h"
#include "ReferenceAngleBondIxn.h"
#include "ReferenceBondForce.h"
#include "ReferenceBrownianDynamics.h"
#include "ReferenceCCMAAlgorithm.h"
#include "ReferenceCMAPTorsionIxn.h"
40
#include "ReferenceConstraints.h"
41
42
#include "ReferenceCustomAngleIxn.h"
#include "ReferenceCustomBondIxn.h"
43
#include "ReferenceCustomCentroidBondIxn.h"
44
#include "ReferenceCustomCompoundBondIxn.h"
45
#include "ReferenceCustomCVForce.h"
46
47
48
49
50
#include "ReferenceCustomDynamics.h"
#include "ReferenceCustomExternalIxn.h"
#include "ReferenceCustomGBIxn.h"
#include "ReferenceCustomHbondIxn.h"
#include "ReferenceCustomNonbondedIxn.h"
51
#include "ReferenceCustomManyParticleIxn.h"
52
#include "ReferenceCustomTorsionIxn.h"
53
#include "ReferenceGayBerneForce.h"
54
55
56
57
58
59
#include "ReferenceHarmonicBondIxn.h"
#include "ReferenceLJCoulomb14.h"
#include "ReferenceLJCoulombIxn.h"
#include "ReferenceMonteCarloBarostat.h"
#include "ReferenceProperDihedralBond.h"
#include "ReferenceRbDihedralBond.h"
60
#include "ReferenceRMSDForce.h"
61
#include "ReferenceStochasticDynamics.h"
62
#include "ReferenceTabulatedFunction.h"
63
64
65
66
#include "ReferenceVariableStochasticDynamics.h"
#include "ReferenceVariableVerletDynamics.h"
#include "ReferenceVerletDynamics.h"
#include "ReferenceVirtualSites.h"
67
#include "openmm/CMMotionRemover.h"
68
#include "openmm/Context.h"
69
#include "openmm/System.h"
70
#include "openmm/internal/AndersenThermostatImpl.h"
71
#include "openmm/internal/ContextImpl.h"
72
#include "openmm/internal/CustomCentroidBondForceImpl.h"
73
#include "openmm/internal/CustomCompoundBondForceImpl.h"
74
#include "openmm/internal/CustomHbondForceImpl.h"
75
#include "openmm/internal/CustomNonbondedForceImpl.h"
76
#include "openmm/internal/CMAPTorsionForceImpl.h"
77
#include "openmm/internal/NonbondedForceImpl.h"
78
#include "openmm/Integrator.h"
79
#include "openmm/OpenMMException.h"
80
#include "SimTKOpenMMUtilities.h"
81
#include "lepton/CustomFunction.h"
82
#include "lepton/Operation.h"
83
84
#include "lepton/Parser.h"
#include "lepton/ParsedExpression.h"
85
#include <cmath>
Peter Eastman's avatar
Peter Eastman committed
86
#include <iostream>
87
#include <limits>
88
89
90
91

using namespace OpenMM;
using namespace std;

peastman's avatar
peastman committed
92
static vector<Vec3>& extractPositions(ContextImpl& context) {
93
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
peastman's avatar
peastman committed
94
    return *((vector<Vec3>*) data->positions);
95
96
}

peastman's avatar
peastman committed
97
static vector<Vec3>& extractVelocities(ContextImpl& context) {
98
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
peastman's avatar
peastman committed
99
    return *((vector<Vec3>*) data->velocities);
100
101
}

peastman's avatar
peastman committed
102
static vector<Vec3>& extractForces(ContextImpl& context) {
103
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
peastman's avatar
peastman committed
104
    return *((vector<Vec3>*) data->forces);
105
106
}

peastman's avatar
peastman committed
107
static Vec3& extractBoxSize(ContextImpl& context) {
108
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
peastman's avatar
peastman committed
109
    return *(Vec3*) data->periodicBoxSize;
110
111
}

peastman's avatar
peastman committed
112
static Vec3* extractBoxVectors(ContextImpl& context) {
113
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
peastman's avatar
peastman committed
114
    return (Vec3*) data->periodicBoxVectors;
115
116
}

117
118
119
120
121
static ReferenceConstraints& extractConstraints(ContextImpl& context) {
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
    return *(ReferenceConstraints*) data->constraints;
}

122
123
124
125
126
static map<string, double>& extractEnergyParameterDerivatives(ContextImpl& context) {
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
    return *((map<string, double>*) data->energyParameterDerivatives);
}

127
128
129
130
131
132
133
/**
 * 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
134
135
    for (auto& child : node.getChildren())
        validateVariables(child, variables);
136
137
}

138
139
140
141
/**
 * Compute the kinetic energy of the system, possibly shifting the velocities in time to account
 * for a leapfrog integrator.
 */
142
static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& masses, double timeShift) {
peastman's avatar
peastman committed
143
144
145
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
146
147
148
149
    int numParticles = context.getSystem().getNumParticles();
    
    // Compute the shifted velocities.
    
peastman's avatar
peastman committed
150
    vector<Vec3> shiftedVel(numParticles);
151
152
153
154
155
    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];
156
    }
157
158
159
    
    // Apply constraints to them.
    
160
161
162
163
    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);
164
165
166
167
168
169
170
    
    // 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]));
171
172
173
    return 0.5*energy;
}

174
void ReferenceCalcForcesAndEnergyKernel::initialize(const System& system) {
175
176
}

177
void ReferenceCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
peastman's avatar
peastman committed
178
    vector<Vec3>& forceData = extractForces(context);
179
180
181
    if (includeForces) {
        int numParticles = context.getSystem().getNumParticles();
        for (int i = 0; i < numParticles; ++i) {
peastman's avatar
peastman committed
182
183
184
            forceData[i][0] = 0.0;
            forceData[i][1] = 0.0;
            forceData[i][2] = 0.0;
185
        }
186
    }
187
188
    else
        savedForces = forceData;
peastman's avatar
peastman committed
189
190
    for (auto& param : context.getParameters())
        extractEnergyParameterDerivatives(context)[param.first] = 0;
191
192
}

193
double ReferenceCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups, bool& valid) {
194
195
    if (!includeForces)
        extractForces(context) = savedForces; // Restore the forces so computing the energy doesn't overwrite the forces with incorrect values.
196
197
    else
        ReferenceVirtualSites::distributeForces(context.getSystem(), extractPositions(context), extractForces(context));
198
199
200
    return 0.0;
}

201
void ReferenceUpdateStateDataKernel::initialize(const System& system) {
202
203
}

204
double ReferenceUpdateStateDataKernel::getTime(const ContextImpl& context) const {
205
206
207
    return data.time;
}

208
void ReferenceUpdateStateDataKernel::setTime(ContextImpl& context, double time) {
209
210
211
    data.time = time;
}

212
213
void ReferenceUpdateStateDataKernel::getPositions(ContextImpl& context, std::vector<Vec3>& positions) {
    int numParticles = context.getSystem().getNumParticles();
peastman's avatar
peastman committed
214
    vector<Vec3>& posData = extractPositions(context);
215
216
217
218
219
220
221
    positions.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        positions[i] = Vec3(posData[i][0], posData[i][1], posData[i][2]);
}

void ReferenceUpdateStateDataKernel::setPositions(ContextImpl& context, const std::vector<Vec3>& positions) {
    int numParticles = context.getSystem().getNumParticles();
peastman's avatar
peastman committed
222
    vector<Vec3>& posData = extractPositions(context);
223
    for (int i = 0; i < numParticles; ++i) {
peastman's avatar
peastman committed
224
225
226
        posData[i][0] = positions[i][0];
        posData[i][1] = positions[i][1];
        posData[i][2] = positions[i][2];
227
228
229
230
231
    }
}

void ReferenceUpdateStateDataKernel::getVelocities(ContextImpl& context, std::vector<Vec3>& velocities) {
    int numParticles = context.getSystem().getNumParticles();
peastman's avatar
peastman committed
232
    vector<Vec3>& velData = extractVelocities(context);
233
234
235
236
237
238
239
    velocities.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        velocities[i] = Vec3(velData[i][0], velData[i][1], velData[i][2]);
}

void ReferenceUpdateStateDataKernel::setVelocities(ContextImpl& context, const std::vector<Vec3>& velocities) {
    int numParticles = context.getSystem().getNumParticles();
peastman's avatar
peastman committed
240
    vector<Vec3>& velData = extractVelocities(context);
241
    for (int i = 0; i < numParticles; ++i) {
peastman's avatar
peastman committed
242
243
244
        velData[i][0] = velocities[i][0];
        velData[i][1] = velocities[i][1];
        velData[i][2] = velocities[i][2];
245
246
247
248
249
    }
}

void ReferenceUpdateStateDataKernel::getForces(ContextImpl& context, std::vector<Vec3>& forces) {
    int numParticles = context.getSystem().getNumParticles();
peastman's avatar
peastman committed
250
    vector<Vec3>& forceData = extractForces(context);
251
252
253
254
255
    forces.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        forces[i] = Vec3(forceData[i][0], forceData[i][1], forceData[i][2]);
}

256
257
258
259
void ReferenceUpdateStateDataKernel::getEnergyParameterDerivatives(ContextImpl& context, map<string, double>& derivs) {
    derivs = extractEnergyParameterDerivatives(context);
}

260
void ReferenceUpdateStateDataKernel::getPeriodicBoxVectors(ContextImpl& context, Vec3& a, Vec3& b, Vec3& c) const {
peastman's avatar
peastman committed
261
    Vec3* vectors = extractBoxVectors(context);
262
263
264
    a = vectors[0];
    b = vectors[1];
    c = vectors[2];
265
266
}

267
void ReferenceUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c) {
peastman's avatar
peastman committed
268
269
270
271
272
    Vec3& box = extractBoxSize(context);
    box[0] = a[0];
    box[1] = b[1];
    box[2] = c[2];
    Vec3* vectors = extractBoxVectors(context);
273
274
275
    vectors[0] = a;
    vectors[1] = b;
    vectors[2] = c;
276
277
}

Peter Eastman's avatar
Peter Eastman committed
278
void ReferenceUpdateStateDataKernel::createCheckpoint(ContextImpl& context, ostream& stream) {
279
    int version = 2;
Peter Eastman's avatar
Peter Eastman committed
280
281
    stream.write((char*) &version, sizeof(int));
    stream.write((char*) &data.time, sizeof(data.time));
peastman's avatar
peastman committed
282
283
284
285
286
287
    vector<Vec3>& posData = extractPositions(context);
    stream.write((char*) &posData[0], sizeof(Vec3)*posData.size());
    vector<Vec3>& velData = extractVelocities(context);
    stream.write((char*) &velData[0], sizeof(Vec3)*velData.size());
    Vec3* vectors = extractBoxVectors(context);
    stream.write((char*) vectors, 3*sizeof(Vec3));
Peter Eastman's avatar
Peter Eastman committed
288
289
290
291
292
293
    SimTKOpenMMUtilities::createCheckpoint(stream);
}

void ReferenceUpdateStateDataKernel::loadCheckpoint(ContextImpl& context, istream& stream) {
    int version;
    stream.read((char*) &version, sizeof(int));
294
    if (version != 2)
Peter Eastman's avatar
Peter Eastman committed
295
296
        throw OpenMMException("Checkpoint was created with a different version of OpenMM");
    stream.read((char*) &data.time, sizeof(data.time));
peastman's avatar
peastman committed
297
298
299
300
301
302
    vector<Vec3>& posData = extractPositions(context);
    stream.read((char*) &posData[0], sizeof(Vec3)*posData.size());
    vector<Vec3>& velData = extractVelocities(context);
    stream.read((char*) &velData[0], sizeof(Vec3)*velData.size());
    Vec3* vectors = extractBoxVectors(context);
    stream.read((char*) vectors, 3*sizeof(Vec3));
Peter Eastman's avatar
Peter Eastman committed
303
304
305
    SimTKOpenMMUtilities::loadCheckpoint(stream);
}

306
307
void ReferenceApplyConstraintsKernel::initialize(const System& system) {
    int numParticles = system.getNumParticles();
308
309
    masses.resize(numParticles);
    inverseMasses.resize(numParticles);
310
    for (int i = 0; i < numParticles; ++i) {
peastman's avatar
peastman committed
311
        masses[i] = system.getParticleMass(i);
312
313
314
315
316
317
318
319
        inverseMasses[i] = 1.0/masses[i];
    }
}

ReferenceApplyConstraintsKernel::~ReferenceApplyConstraintsKernel() {
}

void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
peastman's avatar
peastman committed
320
    vector<Vec3>& positions = extractPositions(context);
321
    extractConstraints(context).apply(positions, positions, inverseMasses, tol);
322
    ReferenceVirtualSites::computePositions(context.getSystem(), positions);
323
324
}

325
void ReferenceApplyConstraintsKernel::applyToVelocities(ContextImpl& context, double tol) {
peastman's avatar
peastman committed
326
327
    vector<Vec3>& positions = extractPositions(context);
    vector<Vec3>& velocities = extractVelocities(context);
328
    extractConstraints(context).applyToVelocities(positions, velocities, inverseMasses, tol);
329
330
}

331
332
333
334
void ReferenceVirtualSitesKernel::initialize(const System& system) {
}

void ReferenceVirtualSitesKernel::computePositions(ContextImpl& context) {
peastman's avatar
peastman committed
335
    vector<Vec3>& positions = extractPositions(context);
336
337
338
    ReferenceVirtualSites::computePositions(context.getSystem(), positions);
}

339
void ReferenceCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
340
    numBonds = force.getNumBonds();
341
342
    bondIndexArray.resize(numBonds, vector<int>(2));
    bondParamArray.resize(numBonds, vector<double>(2));
343
    for (int i = 0; i < numBonds; ++i) {
Peter Eastman's avatar
Peter Eastman committed
344
        int particle1, particle2;
345
        double length, k;
Peter Eastman's avatar
Peter Eastman committed
346
347
348
        force.getBondParameters(i, particle1, particle2, length, k);
        bondIndexArray[i][0] = particle1;
        bondIndexArray[i][1] = particle2;
peastman's avatar
peastman committed
349
350
        bondParamArray[i][0] = length;
        bondParamArray[i][1] = k;
351
    }
352
    usePeriodic = force.usesPeriodicBoundaryConditions();
353
354
}

355
double ReferenceCalcHarmonicBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
356
357
358
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
359
360
    ReferenceBondForce refBondForce;
    ReferenceHarmonicBondIxn harmonicBond;
361
362
    if (usePeriodic)
        harmonicBond.setPeriodic(extractBoxVectors(context));
363
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, includeEnergy ? &energy : NULL, harmonicBond);
364
365
366
    return energy;
}

367
368
369
370
371
372
373
374
375
376
377
378
379
380
void ReferenceCalcHarmonicBondForceKernel::copyParametersToContext(ContextImpl& context, const HarmonicBondForce& force) {
    if (numBonds != force.getNumBonds())
        throw OpenMMException("updateParametersInContext: The number of bonds has changed");

    // Record the values.

    for (int i = 0; i < numBonds; ++i) {
        int particle1, particle2;
        double length, k;
        force.getBondParameters(i, particle1, particle2, length, k);
        if (particle1 != bondIndexArray[i][0] || particle2 != bondIndexArray[i][1])
            throw OpenMMException("updateParametersInContext: The set of particles in a bond has changed");
        bondIndexArray[i][0] = particle1;
        bondIndexArray[i][1] = particle2;
peastman's avatar
peastman committed
381
382
        bondParamArray[i][0] = length;
        bondParamArray[i][1] = k;
383
384
385
    }
}

386
387
388
void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const CustomBondForce& force) {
    numBonds = force.getNumBonds();
    int numParameters = force.getNumPerBondParameters();
389
    usePeriodic = force.usesPeriodicBoundaryConditions();
390
391
392

    // Build the arrays.

393
394
    bondIndexArray.resize(numBonds, vector<int>(2));
    bondParamArray.resize(numBonds, vector<double>(numParameters));
395
    vector<double> params;
396
    for (int i = 0; i < numBonds; ++i) {
397
398
399
400
401
        int particle1, particle2;
        force.getBondParameters(i, particle1, particle2, params);
        bondIndexArray[i][0] = particle1;
        bondIndexArray[i][1] = particle2;
        for (int j = 0; j < numParameters; j++)
peastman's avatar
peastman committed
402
            bondParamArray[i][j] = params[j];
403
404
405
406
407
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
408
409
    energyExpression = expression.createCompiledExpression();
    forceExpression = expression.differentiate("r").createCompiledExpression();
410
411
412
413
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerBondParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
414
415
416
417
418
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
        string param = force.getEnergyParameterDerivativeName(i);
        energyParamDerivNames.push_back(param);
        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
    }
419
420
421
422
423
    set<string> variables;
    variables.insert("r");
    variables.insert(parameterNames.begin(), parameterNames.end());
    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
    validateVariables(expression.getRootNode(), variables);
424
425
}

426
double ReferenceCalcCustomBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
427
428
429
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
430
    map<string, double> globalParameters;
peastman's avatar
peastman committed
431
432
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
433
    ReferenceCustomBondIxn bond(energyExpression, forceExpression, parameterNames, globalParameters, energyParamDerivExpressions);
434
435
    if (usePeriodic)
        bond.setPeriodic(extractBoxVectors(context));
436
437
438
439
440
441
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
    for (int i = 0; i < numBonds; i++)
        bond.calculateBondIxn(bondIndexArray[i], posData, bondParamArray[i], forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
442
443
444
    return energy;
}

445
446
447
448
449
450
451
452
453
454
455
456
457
458
void ReferenceCalcCustomBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomBondForce& force) {
    if (numBonds != force.getNumBonds())
        throw OpenMMException("updateParametersInContext: The number of bonds has changed");

    // Record the values.

    int numParameters = force.getNumPerBondParameters();
    vector<double> params;
    for (int i = 0; i < numBonds; ++i) {
        int particle1, particle2;
        force.getBondParameters(i, particle1, particle2, params);
        if (particle1 != bondIndexArray[i][0] || particle2 != bondIndexArray[i][1])
            throw OpenMMException("updateParametersInContext: The set of particles in a bond has changed");
        for (int j = 0; j < numParameters; j++)
peastman's avatar
peastman committed
459
            bondParamArray[i][j] = params[j];
460
461
462
    }
}

463
464
void ReferenceCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
    numAngles = force.getNumAngles();
465
466
    angleIndexArray.resize(numAngles, vector<int>(3));
    angleParamArray.resize(numAngles, vector<double>(2));
467
    for (int i = 0; i < numAngles; ++i) {
Peter Eastman's avatar
Peter Eastman committed
468
        int particle1, particle2, particle3;
469
        double angle, k;
Peter Eastman's avatar
Peter Eastman committed
470
471
472
473
        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
474
475
        angleParamArray[i][0] = angle;
        angleParamArray[i][1] = k;
476
    }
477
    usePeriodic = force.usesPeriodicBoundaryConditions();
478
479
}

480
double ReferenceCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
481
482
483
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
484
485
    ReferenceBondForce refBondForce;
    ReferenceAngleBondIxn angleBond;
486
487
    if (usePeriodic)
        angleBond.setPeriodic(extractBoxVectors(context));
488
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, includeEnergy ? &energy : NULL, angleBond);
489
490
491
    return energy;
}

492
493
494
495
496
497
498
499
500
501
502
503
void ReferenceCalcHarmonicAngleForceKernel::copyParametersToContext(ContextImpl& context, const HarmonicAngleForce& force) {
    if (numAngles != force.getNumAngles())
        throw OpenMMException("updateParametersInContext: The number of angles has changed");

    // Record the values.

    for (int i = 0; i < numAngles; ++i) {
        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
504
505
        angleParamArray[i][0] = angle;
        angleParamArray[i][1] = k;
506
507
508
    }
}

509
510
511
void ReferenceCalcCustomAngleForceKernel::initialize(const System& system, const CustomAngleForce& force) {
    numAngles = force.getNumAngles();
    int numParameters = force.getNumPerAngleParameters();
512
    usePeriodic = force.usesPeriodicBoundaryConditions();
513
514
515

    // Build the arrays.

516
517
    angleIndexArray.resize(numAngles, vector<int>(3));
    angleParamArray.resize(numAngles, vector<double>(numParameters));
518
    vector<double> params;
519
    for (int i = 0; i < numAngles; ++i) {
520
521
522
523
524
525
        int particle1, particle2, particle3;
        force.getAngleParameters(i, particle1, particle2, particle3, params);
        angleIndexArray[i][0] = particle1;
        angleIndexArray[i][1] = particle2;
        angleIndexArray[i][2] = particle3;
        for (int j = 0; j < numParameters; j++)
peastman's avatar
peastman committed
526
            angleParamArray[i][j] = params[j];
527
528
529
530
531
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
532
533
    energyExpression = expression.createCompiledExpression();
    forceExpression = expression.differentiate("theta").createCompiledExpression();
534
535
536
537
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerAngleParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
538
539
540
541
542
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
        string param = force.getEnergyParameterDerivativeName(i);
        energyParamDerivNames.push_back(param);
        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
    }
543
544
545
546
547
    set<string> variables;
    variables.insert("theta");
    variables.insert(parameterNames.begin(), parameterNames.end());
    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
    validateVariables(expression.getRootNode(), variables);
548
549
}

550
double ReferenceCalcCustomAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
551
552
553
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
554
    map<string, double> globalParameters;
peastman's avatar
peastman committed
555
556
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
557
    ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters, energyParamDerivExpressions);
558
559
    if (usePeriodic)
        customAngle.setPeriodic(extractBoxVectors(context));
560
561
562
563
564
565
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
    for (int i = 0; i < numAngles; i++)
        customAngle.calculateBondIxn(angleIndexArray[i], posData, angleParamArray[i], forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
566
567
568
    return energy;
}

569
570
571
572
573
574
575
576
577
578
579
580
581
582
void ReferenceCalcCustomAngleForceKernel::copyParametersToContext(ContextImpl& context, const CustomAngleForce& force) {
    if (numAngles != force.getNumAngles())
        throw OpenMMException("updateParametersInContext: The number of angles has changed");

    // Record the values.

    int numParameters = force.getNumPerAngleParameters();
    vector<double> params;
    for (int i = 0; i < numAngles; ++i) {
        int particle1, particle2, particle3;
        force.getAngleParameters(i, particle1, particle2, particle3, params);
        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");
        for (int j = 0; j < numParameters; j++)
peastman's avatar
peastman committed
583
            angleParamArray[i][j] = params[j];
584
585
586
    }
}

587
588
void ReferenceCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
    numTorsions = force.getNumTorsions();
589
590
    torsionIndexArray.resize(numTorsions, vector<int>(4));
    torsionParamArray.resize(numTorsions, vector<double>(3));
591
    for (int i = 0; i < numTorsions; ++i) {
Peter Eastman's avatar
Peter Eastman committed
592
        int particle1, particle2, particle3, particle4, periodicity;
593
        double phase, k;
Peter Eastman's avatar
Peter Eastman committed
594
595
596
597
598
        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
599
600
601
        torsionParamArray[i][0] = k;
        torsionParamArray[i][1] = phase;
        torsionParamArray[i][2] = periodicity;
602
    }
603
    usePeriodic = force.usesPeriodicBoundaryConditions();
604
605
}

606
double ReferenceCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
607
608
609
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
610
611
    ReferenceBondForce refBondForce;
    ReferenceProperDihedralBond periodicTorsionBond;
612
613
    if (usePeriodic)
        periodicTorsionBond.setPeriodic(extractBoxVectors(context));
614
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, periodicTorsionBond);
615
616
617
    return energy;
}

618
619
620
621
622
623
624
625
626
627
628
629
void ReferenceCalcPeriodicTorsionForceKernel::copyParametersToContext(ContextImpl& context, const PeriodicTorsionForce& 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, 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
630
631
632
        torsionParamArray[i][0] = k;
        torsionParamArray[i][1] = phase;
        torsionParamArray[i][2] = periodicity;
633
634
635
    }
}

636
637
void ReferenceCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
    numTorsions = force.getNumTorsions();
638
639
    torsionIndexArray.resize(numTorsions, vector<int>(4));
    torsionParamArray.resize(numTorsions, vector<double>(6));
640
    for (int i = 0; i < numTorsions; ++i) {
Peter Eastman's avatar
Peter Eastman committed
641
        int particle1, particle2, particle3, particle4;
642
        double c0, c1, c2, c3, c4, c5;
Peter Eastman's avatar
Peter Eastman committed
643
644
645
646
647
        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
648
649
650
651
652
653
        torsionParamArray[i][0] = c0;
        torsionParamArray[i][1] = c1;
        torsionParamArray[i][2] = c2;
        torsionParamArray[i][3] = c3;
        torsionParamArray[i][4] = c4;
        torsionParamArray[i][5] = c5;
654
    }
655
    usePeriodic = force.usesPeriodicBoundaryConditions();
656
657
}

658
double ReferenceCalcRBTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
659
660
661
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
662
663
    ReferenceBondForce refBondForce;
    ReferenceRbDihedralBond rbTorsionBond;
664
665
    if (usePeriodic)
        rbTorsionBond.setPeriodic(extractBoxVectors(context));
666
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, rbTorsionBond);
667
668
669
    return energy;
}

670
671
672
673
674
675
676
677
678
679
680
681
void ReferenceCalcRBTorsionForceKernel::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
682
683
684
685
686
687
        torsionParamArray[i][0] = c0;
        torsionParamArray[i][1] = c1;
        torsionParamArray[i][2] = c2;
        torsionParamArray[i][3] = c3;
        torsionParamArray[i][4] = c4;
        torsionParamArray[i][5] = c5;
688
689
690
    }
}

691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
void ReferenceCalcCMAPTorsionForceKernel::initialize(const System& system, const CMAPTorsionForce& force) {
    int numMaps = force.getNumMaps();
    int numTorsions = force.getNumTorsions();
    coeff.resize(numMaps);
    vector<double> energy;
    vector<vector<double> > c;
    for (int i = 0; i < numMaps; i++) {
        int size;
        force.getMapParameters(i, size, energy);
        CMAPTorsionForceImpl::calcMapDerivatives(size, energy, c);
        coeff[i].resize(size*size);
        for (int j = 0; j < size*size; j++) {
            coeff[i][j].resize(16);
            for (int k = 0; k < 16; k++)
                coeff[i][j][k] = c[j][k];
        }
    }
    torsionMaps.resize(numTorsions);
    torsionIndices.resize(numTorsions);
    for (int i = 0; i < numTorsions; i++) {
        torsionIndices[i].resize(8);
        force.getTorsionParameters(i, torsionMaps[i], torsionIndices[i][0], torsionIndices[i][1], torsionIndices[i][2],
            torsionIndices[i][3], torsionIndices[i][4], torsionIndices[i][5], torsionIndices[i][6], torsionIndices[i][7]);
    }
715
    usePeriodic = force.usesPeriodicBoundaryConditions();
716
717
}

718
double ReferenceCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
719
720
721
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double totalEnergy = 0;
722
    ReferenceCMAPTorsionIxn torsion(coeff, torsionMaps, torsionIndices);
723
724
    if (usePeriodic)
        torsion.setPeriodic(extractBoxVectors(context));
725
726
727
728
    torsion.calculateIxn(posData, forceData, &totalEnergy);
    return totalEnergy;
}

729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
void ReferenceCalcCMAPTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force) {
    int numMaps = force.getNumMaps();
    int numTorsions = force.getNumTorsions();
    if (coeff.size() != numMaps)
        throw OpenMMException("updateParametersInContext: The number of maps has changed");
    if (torsionMaps.size() != numTorsions)
        throw OpenMMException("updateParametersInContext: The number of CMAP torsions has changed");

    // Update the maps.

    vector<double> energy;
    vector<vector<double> > c;
    for (int i = 0; i < numMaps; i++) {
        int size;
        force.getMapParameters(i, size, energy);
        if (coeff[i].size() != size*size)
            throw OpenMMException("updateParametersInContext: The size of a map has changed");
        CMAPTorsionForceImpl::calcMapDerivatives(size, energy, c);
        for (int j = 0; j < size*size; j++)
            for (int k = 0; k < 16; k++)
                coeff[i][j][k] = c[j][k];
    }

    // Update the indices.

    for (int i = 0; i < numTorsions; i++) {
        int index[8];
        force.getTorsionParameters(i, torsionMaps[i], index[0], index[1], index[2], index[3], index[4], index[5], index[6], index[7]);
        for (int j = 0; j < 8; j++)
            if (index[j] != torsionIndices[i][j])
                throw OpenMMException("updateParametersInContext: The set of particles in a CMAP torsion has changed");
    }
}

763
764
765
void ReferenceCalcCustomTorsionForceKernel::initialize(const System& system, const CustomTorsionForce& force) {
    numTorsions = force.getNumTorsions();
    int numParameters = force.getNumPerTorsionParameters();
766
    usePeriodic = force.usesPeriodicBoundaryConditions();
767
768
769

    // Build the arrays.

770
771
    torsionIndexArray.resize(numTorsions, vector<int>(4));
    torsionParamArray.resize(numTorsions, vector<double>(numParameters));
772
    vector<double> params;
773
    for (int i = 0; i < numTorsions; ++i) {
774
775
776
777
778
779
780
        int particle1, particle2, particle3, particle4;
        force.getTorsionParameters(i, particle1, particle2, particle3, particle4, params);
        torsionIndexArray[i][0] = particle1;
        torsionIndexArray[i][1] = particle2;
        torsionIndexArray[i][2] = particle3;
        torsionIndexArray[i][3] = particle4;
        for (int j = 0; j < numParameters; j++)
peastman's avatar
peastman committed
781
            torsionParamArray[i][j] = params[j];
782
783
784
785
786
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
787
788
    energyExpression = expression.createCompiledExpression();
    forceExpression = expression.differentiate("theta").createCompiledExpression();
789
790
791
792
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerTorsionParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
793
794
795
796
797
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
        string param = force.getEnergyParameterDerivativeName(i);
        energyParamDerivNames.push_back(param);
        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
    }
798
799
800
801
802
    set<string> variables;
    variables.insert("theta");
    variables.insert(parameterNames.begin(), parameterNames.end());
    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
    validateVariables(expression.getRootNode(), variables);
803
804
}

805
double ReferenceCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
806
807
808
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
809
    map<string, double> globalParameters;
peastman's avatar
peastman committed
810
811
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
812
    ReferenceCustomTorsionIxn customTorsion(energyExpression, forceExpression, parameterNames, globalParameters, energyParamDerivExpressions);
813
814
    if (usePeriodic)
        customTorsion.setPeriodic(extractBoxVectors(context));
815
816
817
818
819
820
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
    for (int i = 0; i < numTorsions; i++)
        customTorsion.calculateBondIxn(torsionIndexArray[i], posData, torsionParamArray[i], forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
821
822
823
    return energy;
}

824
825
826
827
828
829
830
831
832
833
834
835
836
837
void ReferenceCalcCustomTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CustomTorsionForce& force) {
    if (numTorsions != force.getNumTorsions())
        throw OpenMMException("updateParametersInContext: The number of torsions has changed");

    // Record the values.

    int numParameters = force.getNumPerTorsionParameters();
    vector<double> params;
    for (int i = 0; i < numTorsions; ++i) {
        int particle1, particle2, particle3, particle4;
        force.getTorsionParameters(i, particle1, particle2, particle3, particle4, params);
        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");
        for (int j = 0; j < numParameters; j++)
peastman's avatar
peastman committed
838
            torsionParamArray[i][j] = params[j];
839
840
841
    }
}

842
843
844
845
846
ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
    if (neighborList != NULL)
        delete neighborList;
}

847
848
849
850
void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {

    // Identify which exceptions are 1-4 interactions.

Peter Eastman's avatar
Peter Eastman committed
851
    numParticles = force.getNumParticles();
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
    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);
        if (chargeProd != 0.0 || epsilon != 0.0)
            nb14s.push_back(i);
    }

    // Build the arrays.

    num14 = nb14s.size();
867
868
869
    bonded14IndexArray.resize(num14, vector<int>(2));
    bonded14ParamArray.resize(num14, vector<double>(3));
    particleParamArray.resize(numParticles, vector<double>(3));
870
871
872
873
    baseParticleParams.resize(numParticles);
    baseExceptionParams.resize(num14);
    for (int i = 0; i < numParticles; ++i)
       force.getParticleParameters(i, baseParticleParams[i][0], baseParticleParams[i][1], baseParticleParams[i][2]);
874
    this->exclusions = exclusions;
875
    for (int i = 0; i < num14; ++i) {
Peter Eastman's avatar
Peter Eastman committed
876
        int particle1, particle2;
877
        force.getExceptionParameters(nb14s[i], particle1, particle2, baseExceptionParams[i][0], baseExceptionParams[i][1], baseExceptionParams[i][2]);
Peter Eastman's avatar
Peter Eastman committed
878
879
        bonded14IndexArray[i][0] = particle1;
        bonded14IndexArray[i][1] = particle2;
880
881
882
883
884
885
886
887
888
889
890
891
892
893
    }
    for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
        string param;
        int particle;
        double charge, sigma, epsilon;
        force.getParticleParameterOffset(i, param, particle, charge, sigma, epsilon);
        particleParamOffsets[make_pair(param, particle)] = {charge, sigma, epsilon};
    }
    for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
        string param;
        int exception;
        double charge, sigma, epsilon;
        force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
        exceptionParamOffsets[make_pair(param, exception)] = {charge, sigma, epsilon};
894
    }
895
    nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
peastman's avatar
peastman committed
896
    nonbondedCutoff = force.getCutoffDistance();
897
    if (nonbondedMethod == NoCutoff) {
898
        neighborList = NULL;
899
900
901
        useSwitchingFunction = false;
    }
    else {
902
        neighborList = new NeighborList();
903
904
905
        useSwitchingFunction = force.getUseSwitchingFunction();
        switchingDistance = force.getSwitchingDistance();
    }
906
907
908
    if (nonbondedMethod == Ewald) {
        double alpha;
        NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmax[0], kmax[1], kmax[2]);
peastman's avatar
peastman committed
909
        ewaldAlpha = alpha;
910
911
912
    }
    else if (nonbondedMethod == PME) {
        double alpha;
913
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2], false);
peastman's avatar
peastman committed
914
        ewaldAlpha = alpha;
915
    }
916
917
918
    else if (nonbondedMethod == LJPME) {
        double alpha;
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2], false);
peastman's avatar
peastman committed
919
        ewaldAlpha = alpha;
920
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, dispersionGridSize[0], dispersionGridSize[1], dispersionGridSize[2], true);
peastman's avatar
peastman committed
921
        ewaldDispersionAlpha = alpha;
922
923
        useSwitchingFunction = false;
    }
peastman's avatar
peastman committed
924
    rfDielectric = force.getReactionFieldDielectric();
925
926
927
928
    if (force.getUseDispersionCorrection())
        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
    else
        dispersionCoefficient = 0.0;
929
930
}

931
double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
932
    computeParameters(context);
peastman's avatar
peastman committed
933
934
935
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
936
    ReferenceLJCoulombIxn clj;
937
    bool periodic = (nonbondedMethod == CutoffPeriodic);
938
    bool ewald  = (nonbondedMethod == Ewald);
939
    bool pme  = (nonbondedMethod == PME);
940
    bool ljpme = (nonbondedMethod == LJPME);
941
    if (nonbondedMethod != NoCutoff) {
942
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic || ewald || pme || ljpme, nonbondedCutoff, 0.0);
943
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
944
    }
945
    if (periodic || ewald || pme || ljpme) {
peastman's avatar
peastman committed
946
        Vec3* boxVectors = extractBoxVectors(context);
947
        double minAllowedSize = 1.999999*nonbondedCutoff;
948
        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
949
            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
950
        clj.setPeriodic(boxVectors);
951
    }
952
953
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
954
    if (pme)
955
        clj.setUsePME(ewaldAlpha, gridSize);
956
957
958
959
    if (ljpme){
        clj.setUsePME(ewaldAlpha, gridSize);
        clj.setUseLJPME(ewaldDispersionAlpha, dispersionGridSize);
    }
960
961
    if (useSwitchingFunction)
        clj.setUseSwitchingFunction(switchingDistance);
962
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusions, forceData, includeEnergy ? &energy : NULL, includeDirect, includeReciprocal);
963
964
965
966
    if (includeDirect) {
        ReferenceBondForce refBondForce;
        ReferenceLJCoulomb14 nonbonded14;
        refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
967
        if (periodic || ewald || pme) {
peastman's avatar
peastman committed
968
            Vec3* boxVectors = extractBoxVectors(context);
969
            energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
970
        }
971
    }
972
973
974
    return energy;
}

975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
void ReferenceCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const NonbondedForce& force) {
    if (force.getNumParticles() != numParticles)
        throw OpenMMException("updateParametersInContext: The number of particles has changed");
    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);
        if (chargeProd != 0.0 || epsilon != 0.0)
            nb14s.push_back(i);
    }
    if (nb14s.size() != num14)
        throw OpenMMException("updateParametersInContext: The number of non-excluded exceptions has changed");

    // Record the values.

991
992
    for (int i = 0; i < numParticles; ++i)
        force.getParticleParameters(i, baseParticleParams[i][0], baseParticleParams[i][1], baseParticleParams[i][2]);
993
994
    for (int i = 0; i < num14; ++i) {
        int particle1, particle2;
995
        force.getExceptionParameters(nb14s[i], particle1, particle2, baseExceptionParams[i][0], baseExceptionParams[i][1], baseExceptionParams[i][2]);
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
        bonded14IndexArray[i][0] = particle1;
        bonded14IndexArray[i][1] = particle2;
    }
    
    // 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);
}

1007
void ReferenceCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
1008
1009
    if (nonbondedMethod != PME && nonbondedMethod != LJPME)
        throw OpenMMException("getPMEParametersInContext: This Context is not using PME or LJPME");
1010
1011
1012
1013
1014
1015
    alpha = ewaldAlpha;
    nx = gridSize[0];
    ny = gridSize[1];
    nz = gridSize[2];
}

1016
void ReferenceCalcNonbondedForceKernel::getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
1017
1018
    if (nonbondedMethod != LJPME)
        throw OpenMMException("getPMEParametersInContext: This Context is not using LJPME");
1019
1020
1021
1022
    alpha = ewaldDispersionAlpha;
    nx = dispersionGridSize[0];
    ny = dispersionGridSize[1];
    nz = dispersionGridSize[2];
1023
1024
}

1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
void ReferenceCalcNonbondedForceKernel::computeParameters(ContextImpl& context) {
    // Compute particle parameters.

    vector<double> charges(numParticles), sigmas(numParticles), epsilons(numParticles);
    for (int i = 0; i < numParticles; i++) {
        charges[i] = baseParticleParams[i][0];
        sigmas[i] = baseParticleParams[i][1];
        epsilons[i] = baseParticleParams[i][2];
    }
    for (auto& offset : particleParamOffsets) {
        double value = context.getParameter(offset.first.first);
        int index = offset.first.second;
        charges[index] += value*offset.second[0];
        sigmas[index] += value*offset.second[1];
        epsilons[index] += value*offset.second[2];
    }
    for (int i = 0; i < numParticles; i++) {
        particleParamArray[i][0] = 0.5*sigmas[i];
        particleParamArray[i][1] = 2.0*sqrt(epsilons[i]);
        particleParamArray[i][2] = charges[i];
    }

    // Compute exception parameters.

    charges.resize(num14);
    sigmas.resize(num14);
    epsilons.resize(num14);
    for (int i = 0; i < num14; i++) {
        charges[i] = baseExceptionParams[i][0];
        sigmas[i] = baseExceptionParams[i][1];
        epsilons[i] = baseExceptionParams[i][2];
    }
    for (auto& offset : exceptionParamOffsets) {
        double value = context.getParameter(offset.first.first);
        int index = offset.first.second;
        charges[index] += value*offset.second[0];
        sigmas[index] += value*offset.second[1];
        epsilons[index] += value*offset.second[2];
    }
    for (int i = 0; i < num14; i++) {
        bonded14ParamArray[i][0] = sigmas[i];
        bonded14ParamArray[i][1] = 4.0*epsilons[i];
        bonded14ParamArray[i][2] = charges[i];
    }
}

1071
1072
1073
ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
    if (neighborList != NULL)
        delete neighborList;
1074
1075
    if (forceCopy != NULL)
        delete forceCopy;
1076
1077
1078
1079
}

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

1080
    // Record the exclusions.
1081
1082
1083

    numParticles = force.getNumParticles();
    exclusions.resize(numParticles);
1084
    for (int i = 0; i < force.getNumExclusions(); i++) {
1085
        int particle1, particle2;
1086
        force.getExclusionParticles(i, particle1, particle2);
1087
1088
1089
1090
1091
1092
        exclusions[particle1].insert(particle2);
        exclusions[particle2].insert(particle1);
    }

    // Build the arrays.

1093
    int numParameters = force.getNumPerParticleParameters();
1094
1095
1096
    particleParamArray.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        force.getParticleParameters(i, particleParamArray[i]);
1097
    nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
peastman's avatar
peastman committed
1098
    nonbondedCutoff = force.getCutoffDistance();
1099
    if (nonbondedMethod == NoCutoff) {
1100
        neighborList = NULL;
1101
1102
1103
        useSwitchingFunction = false;
    }
    else {
1104
        neighborList = new NeighborList();
1105
1106
1107
        useSwitchingFunction = force.getUseSwitchingFunction();
        switchingDistance = force.getSwitchingDistance();
    }
1108

1109
1110
1111
    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
1112
    for (int i = 0; i < force.getNumFunctions(); i++)
1113
        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
1114

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

1117
    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
1118
1119
    energyExpression = expression.createCompiledExpression();
    forceExpression = expression.differentiate("r").createCompiledExpression();
1120
1121
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerParticleParameterName(i));
1122
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
1123
        globalParameterNames.push_back(force.getGlobalParameterName(i));
1124
1125
        globalParamValues[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
    }
1126
1127
1128
1129
1130
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
        string param = force.getEnergyParameterDerivativeName(i);
        energyParamDerivNames.push_back(param);
        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
    }
1131
1132
1133
1134
1135
1136
1137
1138
    set<string> variables;
    variables.insert("r");
    for (int i = 0; i < numParameters; i++) {
        variables.insert(parameterNames[i]+"1");
        variables.insert(parameterNames[i]+"2");
    }
    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
    validateVariables(expression.getRootNode(), variables);
1139
1140
1141

    // Delete the custom functions.

peastman's avatar
peastman committed
1142
1143
    for (auto& function : functions)
        delete function.second;
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
    
    // 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;
    }
1155
1156
1157
1158
1159
1160
1161
1162
    
    // 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));
    }
1163
1164
}

1165
double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1166
1167
1168
1169
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    Vec3* boxVectors = extractBoxVectors(context);
    double energy = 0;
1170
    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, energyParamDerivExpressions);
1171
1172
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
1173
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic, nonbondedCutoff, 0.0);
1174
1175
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
1176
1177
    if (periodic) {
        double minAllowedSize = 2*nonbondedCutoff;
1178
        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
1179
            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
1180
        ixn.setPeriodic(boxVectors);
1181
    }
1182
1183
    if (interactionGroups.size() > 0)
        ixn.setInteractionGroups(interactionGroups);
1184
    bool globalParamsChanged = false;
peastman's avatar
peastman committed
1185
1186
1187
    for (auto& name : globalParameterNames) {
        double value = context.getParameter(name);
        if (globalParamValues[name] != value)
1188
            globalParamsChanged = true;
peastman's avatar
peastman committed
1189
        globalParamValues[name] = value;
1190
    }
1191
1192
    if (useSwitchingFunction)
        ixn.setUseSwitchingFunction(switchingDistance);
1193
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
1194
    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusions, globalParamValues, forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
1195
1196
1197
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
1198
1199
1200
1201
    
    // Add in the long range correction.
    
    if (!hasInitializedLongRangeCorrection || (globalParamsChanged && forceCopy != NULL)) {
1202
        CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, context.getOwner(), longRangeCoefficient, longRangeCoefficientDerivs);
1203
1204
        hasInitializedLongRangeCorrection = true;
    }
1205
1206
1207
1208
    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;
1209
1210
1211
    return energy;
}

1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
void ReferenceCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& 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
1224
            particleParamArray[i][j] = parameters[j];
1225
    }
1226
1227
1228
1229
    
    // If necessary, recompute the long range correction.
    
    if (forceCopy != NULL) {
1230
        CustomNonbondedForceImpl::calcLongRangeCorrection(force, context.getOwner(), longRangeCoefficient, longRangeCoefficientDerivs);
1231
1232
1233
        hasInitializedLongRangeCorrection = true;
        *forceCopy = force;
    }
1234
1235
}

1236
ReferenceCalcGBSAOBCForceKernel::~ReferenceCalcGBSAOBCForceKernel() {
1237
    if (obc) {
Peter Eastman's avatar
Peter Eastman committed
1238
        delete obc->getObcParameters();
1239
1240
1241
1242
        delete obc;
    }
}

1243
void ReferenceCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
Peter Eastman's avatar
Peter Eastman committed
1244
1245
    int numParticles = system.getNumParticles();
    charges.resize(numParticles);
peastman's avatar
peastman committed
1246
1247
    vector<double> atomicRadii(numParticles);
    vector<double> scaleFactors(numParticles);
Peter Eastman's avatar
Peter Eastman committed
1248
    for (int i = 0; i < numParticles; ++i) {
1249
        double charge, radius, scalingFactor;
Peter Eastman's avatar
Peter Eastman committed
1250
        force.getParticleParameters(i, charge, radius, scalingFactor);
peastman's avatar
peastman committed
1251
1252
1253
        charges[i] = charge;
        atomicRadii[i] = radius;
        scaleFactors[i] = scalingFactor;
1254
    }
1255
    ObcParameters* obcParameters = new ObcParameters(numParticles, ObcParameters::ObcTypeII);
1256
    obcParameters->setAtomicRadii(atomicRadii);
1257
    obcParameters->setScaledRadiusFactors(scaleFactors);
peastman's avatar
peastman committed
1258
1259
    obcParameters->setSolventDielectric(force.getSolventDielectric());
    obcParameters->setSoluteDielectric(force.getSoluteDielectric());
1260
    obcParameters->setPi4Asolv(4*M_PI*force.getSurfaceAreaEnergy());
1261
    if (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff)
peastman's avatar
peastman committed
1262
        obcParameters->setUseCutoff(force.getCutoffDistance());
1263
    isPeriodic = (force.getNonbondedMethod() == GBSAOBCForce::CutoffPeriodic);
1264
    obc = new ReferenceObc(obcParameters);
1265
    obc->setIncludeAceApproximation(true);
1266
1267
}

1268
double ReferenceCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1269
1270
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
1271
    if (isPeriodic)
1272
        obc->getObcParameters()->setPeriodic(extractBoxVectors(context));
Mark Friedrichs's avatar
Mark Friedrichs committed
1273
    return obc->computeBornEnergyForces(posData, charges, forceData);
1274
1275
}

1276
1277
1278
1279
1280
1281
1282
1283
void ReferenceCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force) {
    int numParticles = force.getNumParticles();
    ObcParameters* obcParameters = obc->getObcParameters();
    if (numParticles != obcParameters->getAtomicRadii().size())
        throw OpenMMException("updateParametersInContext: The number of particles has changed");

    // Record the values.

peastman's avatar
peastman committed
1284
1285
    vector<double> atomicRadii(numParticles);
    vector<double> scaleFactors(numParticles);
1286
1287
1288
    for (int i = 0; i < numParticles; ++i) {
        double charge, radius, scalingFactor;
        force.getParticleParameters(i, charge, radius, scalingFactor);
peastman's avatar
peastman committed
1289
1290
1291
        charges[i] = charge;
        atomicRadii[i] = radius;
        scaleFactors[i] = scalingFactor;
1292
1293
1294
1295
1296
    }
    obcParameters->setAtomicRadii(atomicRadii);
    obcParameters->setScaledRadiusFactors(scaleFactors);
}

1297
1298
1299
1300
1301
1302
ReferenceCalcCustomGBForceKernel::~ReferenceCalcCustomGBForceKernel() {
    if (neighborList != NULL)
        delete neighborList;
}

void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
    if (force.getNumComputedValues() > 0) {
        string name, expression;
        CustomGBForce::ComputationType type;
        force.getComputedValueParameters(0, name, expression, type);
        if (type == CustomGBForce::SingleParticle)
            throw OpenMMException("ReferencePlatform 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("ReferencePlatform requires that a CustomGBForce only have one computed value of type ParticlePair or ParticlePairNoExclusions.");
        }
    }
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329

    // 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.

    int numPerParticleParameters = force.getNumPerParticleParameters();
1330
1331
1332
    particleParamArray.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        force.getParticleParameters(i, particleParamArray[i]);
1333
1334
1335
1336
1337
    for (int i = 0; i < numPerParticleParameters; i++)
        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
1338
    nonbondedCutoff = force.getCutoffDistance();
1339
1340
1341
1342
1343
1344
1345
1346
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();

    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
1347
    for (int i = 0; i < force.getNumFunctions(); i++)
1348
        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
1349
1350
1351

    // Parse the expressions for computed values.

1352
    valueDerivExpressions.resize(force.getNumComputedValues());
1353
    valueGradientExpressions.resize(force.getNumComputedValues());
1354
    valueParamDerivExpressions.resize(force.getNumComputedValues());
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
    set<string> particleVariables, pairVariables;
    pairVariables.insert("r");
    particleVariables.insert("x");
    particleVariables.insert("y");
    particleVariables.insert("z");
    for (int i = 0; i < numPerParticleParameters; i++) {
        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());
1367
1368
1369
1370
1371
    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();
1372
        valueExpressions.push_back(ex.createCompiledExpression());
1373
1374
        valueTypes.push_back(type);
        valueNames.push_back(name);
1375
        if (i == 0) {
1376
            valueDerivExpressions[i].push_back(ex.differentiate("r").createCompiledExpression());
1377
1378
            validateVariables(ex.getRootNode(), pairVariables);
        }
1379
        else {
1380
1381
1382
            valueGradientExpressions[i].push_back(ex.differentiate("x").createCompiledExpression());
            valueGradientExpressions[i].push_back(ex.differentiate("y").createCompiledExpression());
            valueGradientExpressions[i].push_back(ex.differentiate("z").createCompiledExpression());
1383
            for (int j = 0; j < i; j++)
1384
                valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).createCompiledExpression());
1385
            validateVariables(ex.getRootNode(), particleVariables);
1386
        }
1387
1388
1389
1390
1391
        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());
        }
1392
1393
1394
        particleVariables.insert(name);
        pairVariables.insert(name+"1");
        pairVariables.insert(name+"2");
1395
1396
    }

1397
    // Parse the expressions for energy terms.
1398
1399

    energyDerivExpressions.resize(force.getNumEnergyTerms());
1400
    energyGradientExpressions.resize(force.getNumEnergyTerms());
1401
    energyParamDerivExpressions.resize(force.getNumEnergyTerms());
1402
1403
1404
1405
1406
    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();
1407
        energyExpressions.push_back(ex.createCompiledExpression());
1408
1409
        energyTypes.push_back(type);
        if (type != CustomGBForce::SingleParticle)
1410
            energyDerivExpressions[i].push_back(ex.differentiate("r").createCompiledExpression());
1411
        for (int j = 0; j < force.getNumComputedValues(); j++) {
1412
            if (type == CustomGBForce::SingleParticle) {
1413
1414
1415
1416
                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());
1417
                validateVariables(ex.getRootNode(), particleVariables);
1418
            }
1419
            else {
1420
1421
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").createCompiledExpression());
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").createCompiledExpression());
1422
                validateVariables(ex.getRootNode(), pairVariables);
1423
1424
            }
        }
1425
1426
        for (int j = 0; j < force.getNumEnergyParameterDerivatives(); j++)
            energyParamDerivExpressions[i].push_back(ex.differentiate(force.getEnergyParameterDerivativeName(j)).createCompiledExpression());
1427
1428
1429
1430
    }

    // Delete the custom functions.

peastman's avatar
peastman committed
1431
1432
    for (auto& function : functions)
        delete function.second;
1433
1434
}

1435
double ReferenceCalcCustomGBForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1436
1437
1438
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
1439
1440
    ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueParamDerivExpressions, valueNames, valueTypes,
        energyExpressions, energyDerivExpressions, energyGradientExpressions, energyParamDerivExpressions, energyTypes, particleParameterNames);
1441
    bool periodic = (nonbondedMethod == CutoffPeriodic);
1442
    if (periodic)
1443
        ixn.setPeriodic(extractBoxVectors(context));
1444
    if (nonbondedMethod != NoCutoff) {
Peter Eastman's avatar
Peter Eastman committed
1445
1446
        vector<set<int> > empty(context.getSystem().getNumParticles()); // Don't omit exclusions from the neighbor list
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, empty, extractBoxVectors(context), periodic, nonbondedCutoff, 0.0);
1447
1448
1449
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1450
1451
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1452
1453
1454
1455
1456
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
    ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
1457
1458
1459
    return energy;
}

1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
void ReferenceCalcCustomGBForceKernel::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
1472
            particleParamArray[i][j] = parameters[j];
1473
1474
1475
    }
}

peastman's avatar
peastman committed
1476
ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::PeriodicDistanceFunction(Vec3** boxVectorHandle) : boxVectorHandle(boxVectorHandle) {
1477
1478
1479
1480
1481
1482
1483
}

int ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::getNumArguments() const {
    return 6;
}

double ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::evaluate(const double* arguments) const {
peastman's avatar
peastman committed
1484
1485
    Vec3* boxVectors = *boxVectorHandle;
    Vec3 delta = Vec3(arguments[0], arguments[1], arguments[2])-Vec3(arguments[3], arguments[4], arguments[5]);
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
    delta -= boxVectors[2]*floor(delta[2]/boxVectors[2][2]+0.5);
    delta -= boxVectors[1]*floor(delta[1]/boxVectors[1][1]+0.5);
    delta -= boxVectors[0]*floor(delta[0]/boxVectors[0][0]+0.5);
    return sqrt(delta.dot(delta));
}

double ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
    int argIndex = -1;
    for (int i = 0; i < 6; i++) {
        if (derivOrder[i] > 0) {
            if (derivOrder[i] > 1 || argIndex != -1)
                throw OpenMMException("Unsupported derivative of periodicdistance"); // Should be impossible for this to happen.
            argIndex = i;
        }
    }
peastman's avatar
peastman committed
1501
1502
    Vec3* boxVectors = *boxVectorHandle;
    Vec3 delta = Vec3(arguments[0], arguments[1], arguments[2])-Vec3(arguments[3], arguments[4], arguments[5]);
1503
1504
1505
1506
    delta -= boxVectors[2]*floor(delta[2]/boxVectors[2][2]+0.5);
    delta -= boxVectors[1]*floor(delta[1]/boxVectors[1][1]+0.5);
    delta -= boxVectors[0]*floor(delta[0]/boxVectors[0][0]+0.5);
    double r = sqrt(delta.dot(delta));
1507
1508
    if (r == 0)
        return 0.0;    
1509
1510
1511
1512
1513
1514
1515
1516
1517
    if (argIndex < 3)
        return delta[argIndex]/r;
    return -delta[argIndex-3]/r;
}

Lepton::CustomFunction* ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::clone() const {
    return new PeriodicDistanceFunction(boxVectorHandle);
}

1518
1519
1520
1521
1522
1523
1524
void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, const CustomExternalForce& force) {
    numParticles = force.getNumParticles();
    int numParameters = force.getNumPerParticleParameters();

    // Build the arrays.

    particles.resize(numParticles);
1525
1526
1527
    particleParamArray.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        force.getParticleParameters(i, particles[i], particleParamArray[i]);
1528
1529
1530

    // Parse the expression used to calculate the force.

1531
1532
1533
1534
    map<string, Lepton::CustomFunction*> functions;
    PeriodicDistanceFunction periodicDistance(&boxVectors);
    functions["periodicdistance"] = &periodicDistance;
    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
1535
1536
1537
1538
    energyExpression = expression.createCompiledExpression();
    forceExpressionX = expression.differentiate("x").createCompiledExpression();
    forceExpressionY = expression.differentiate("y").createCompiledExpression();
    forceExpressionZ = expression.differentiate("z").createCompiledExpression();
1539
1540
1541
1542
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerParticleParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
1543
1544
1545
1546
1547
1548
1549
    set<string> variables;
    variables.insert("x");
    variables.insert("y");
    variables.insert("z");
    variables.insert(parameterNames.begin(), parameterNames.end());
    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
    validateVariables(expression.getRootNode(), variables);
1550
1551
}

1552
double ReferenceCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1553
1554
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
1555
    boxVectors = extractBoxVectors(context);
peastman's avatar
peastman committed
1556
    double energy = 0;
1557
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1558
1559
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1560
1561
    ReferenceCustomExternalIxn force(energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ, parameterNames, globalParameters);
    for (int i = 0; i < numParticles; ++i)
1562
        force.calculateForce(particles[i], posData, particleParamArray[i], forceData, includeEnergy ? &energy : NULL);
1563
1564
1565
    return energy;
}

1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
void ReferenceCalcCustomExternalForceKernel::copyParametersToContext(ContextImpl& context, const CustomExternalForce& 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) {
        int particle;
        vector<double> parameters;
        force.getParticleParameters(i, particle, parameters);
        if (particle != particles[i])
            throw OpenMMException("updateParametersInContext: A particle index has changed");
        for (int j = 0; j < numParameters; j++)
peastman's avatar
peastman committed
1581
            particleParamArray[i][j] = parameters[j];
1582
1583
1584
    }
}

1585
ReferenceCalcCustomHbondForceKernel::~ReferenceCalcCustomHbondForceKernel() {
1586
1587
    if (ixn != NULL)
        delete ixn;
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
}

void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const CustomHbondForce& force) {

    // Record the exclusions.

    numDonors = force.getNumDonors();
    numAcceptors = force.getNumAcceptors();
    numParticles = system.getNumParticles();
    exclusions.resize(numDonors);
    for (int i = 0; i < force.getNumExclusions(); i++) {
        int donor, acceptor;
        force.getExclusionParticles(i, donor, acceptor);
        exclusions[donor].insert(acceptor);
    }

    // Build the arrays.

1606
    vector<vector<int> > donorParticles(numDonors);
1607
    int numDonorParameters = force.getNumPerDonorParameters();
1608
    donorParamArray.resize(numDonors);
1609
    for (int i = 0; i < numDonors; ++i) {
1610
        int d1, d2, d3;
1611
        force.getDonorParameters(i, d1, d2, d3, donorParamArray[i]);
1612
1613
1614
        donorParticles[i].push_back(d1);
        donorParticles[i].push_back(d2);
        donorParticles[i].push_back(d3);
1615
    }
1616
    vector<vector<int> > acceptorParticles(numAcceptors);
1617
    int numAcceptorParameters = force.getNumPerAcceptorParameters();
1618
    acceptorParamArray.resize(numAcceptors);
1619
    for (int i = 0; i < numAcceptors; ++i) {
1620
        int a1, a2, a3;
1621
        force.getAcceptorParameters(i, a1, a2, a3, acceptorParamArray[i]);
1622
1623
1624
        acceptorParticles[i].push_back(a1);
        acceptorParticles[i].push_back(a2);
        acceptorParticles[i].push_back(a3);
1625
    }
1626
    NonbondedMethod nonbondedMethod = CalcCustomHbondForceKernel::NonbondedMethod(force.getNonbondedMethod());
peastman's avatar
peastman committed
1627
    nonbondedCutoff = force.getCutoffDistance();
1628
1629
1630
1631

    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
1632
    for (int i = 0; i < force.getNumFunctions(); i++)
1633
        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
1634

1635
    // Parse the expression and create the object used to calculate the interaction.
1636

1637
1638
1639
    map<string, vector<int> > distances;
    map<string, vector<int> > angles;
    map<string, vector<int> > dihedrals;
1640
    Lepton::ParsedExpression energyExpression = CustomHbondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
1641
1642
    vector<string> donorParameterNames;
    vector<string> acceptorParameterNames;
1643
1644
1645
1646
1647
1648
    for (int i = 0; i < numDonorParameters; i++)
        donorParameterNames.push_back(force.getPerDonorParameterName(i));
    for (int i = 0; i < numAcceptorParameters; i++)
        acceptorParameterNames.push_back(force.getPerAcceptorParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
1649
    ixn = new ReferenceCustomHbondIxn(donorParticles, acceptorParticles, energyExpression, donorParameterNames, acceptorParameterNames, distances, angles, dihedrals);
1650
    isPeriodic = (nonbondedMethod == CutoffPeriodic);
1651
1652
    if (nonbondedMethod != NoCutoff)
        ixn->setUseCutoff(nonbondedCutoff);
1653
1654
1655

    // Delete the custom functions.

peastman's avatar
peastman committed
1656
1657
    for (auto& function : functions)
        delete function.second;
1658
1659
}

1660
double ReferenceCalcCustomHbondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1661
1662
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
1663
    if (isPeriodic)
1664
        ixn->setPeriodic(extractBoxVectors(context));
peastman's avatar
peastman committed
1665
    double energy = 0;
1666
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1667
1668
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1669
    ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusions, globalParameters, forceData, includeEnergy ? &energy : NULL);
1670
1671
1672
    return energy;
}

1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
void ReferenceCalcCustomHbondForceKernel::copyParametersToContext(ContextImpl& context, const CustomHbondForce& force) {
    if (numDonors != force.getNumDonors())
        throw OpenMMException("updateParametersInContext: The number of donors has changed");
    if (numAcceptors != force.getNumAcceptors())
        throw OpenMMException("updateParametersInContext: The number of acceptors has changed");

    // Record the values.

    vector<double> parameters;
    int numDonorParameters = force.getNumPerDonorParameters();
    const vector<vector<int> >& donorAtoms = ixn->getDonorAtoms();
    for (int i = 0; i < numDonors; ++i) {
        int d1, d2, d3;
        force.getDonorParameters(i, d1, d2, d3, parameters);
        if (d1 != donorAtoms[i][0] || d2 != donorAtoms[i][1] || d3 != donorAtoms[i][2])
            throw OpenMMException("updateParametersInContext: The set of particles in a donor group has changed");
        for (int j = 0; j < numDonorParameters; j++)
peastman's avatar
peastman committed
1690
            donorParamArray[i][j] = parameters[j];
1691
1692
1693
1694
1695
1696
1697
1698
1699
    }
    int numAcceptorParameters = force.getNumPerAcceptorParameters();
    const vector<vector<int> >& acceptorAtoms = ixn->getAcceptorAtoms();
    for (int i = 0; i < numAcceptors; ++i) {
        int a1, a2, a3;
        force.getAcceptorParameters(i, a1, a2, a3, parameters);
        if (a1 != acceptorAtoms[i][0] || a2 != acceptorAtoms[i][1] || a3 != acceptorAtoms[i][2])
            throw OpenMMException("updateParametersInContext: The set of particles in an acceptor group has changed");
        for (int j = 0; j < numAcceptorParameters; j++)
peastman's avatar
peastman committed
1700
            acceptorParamArray[i][j] = parameters[j];
1701
1702
1703
    }
}

1704
1705
1706
1707
1708
1709
ReferenceCalcCustomCentroidBondForceKernel::~ReferenceCalcCustomCentroidBondForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

void ReferenceCalcCustomCentroidBondForceKernel::initialize(const System& system, const CustomCentroidBondForce& force) {
1710
    usePeriodic = force.usesPeriodicBoundaryConditions();
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723

    // Build the arrays.

    int numGroups = force.getNumGroups();
    vector<vector<int> > groupAtoms(numGroups);
    vector<double> ignored;
    for (int i = 0; i < numGroups; i++)
        force.getGroupParameters(i, groupAtoms[i], ignored);
    vector<vector<double> > normalizedWeights;
    CustomCentroidBondForceImpl::computeNormalizedWeights(force, system, normalizedWeights);
    numBonds = force.getNumBonds();
    vector<vector<int> > bondGroups(numBonds);
    int numBondParameters = force.getNumPerBondParameters();
1724
1725
1726
    bondParamArray.resize(numBonds);
    for (int i = 0; i < numBonds; ++i)
        force.getBondParameters(i, bondGroups[i], bondParamArray[i]);
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744

    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
    for (int i = 0; i < force.getNumFunctions(); i++)
        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));

    // Parse the expression and create the object used to calculate the interaction.

    map<string, vector<int> > distances;
    map<string, vector<int> > angles;
    map<string, vector<int> > dihedrals;
    Lepton::ParsedExpression energyExpression = CustomCentroidBondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
    vector<string> bondParameterNames;
    for (int i = 0; i < numBondParameters; i++)
        bondParameterNames.push_back(force.getPerBondParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
1745
1746
1747
1748
1749
1750
1751
    vector<Lepton::CompiledExpression> energyParamDerivExpressions;
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
        string param = force.getEnergyParameterDerivativeName(i);
        energyParamDerivNames.push_back(param);
        energyParamDerivExpressions.push_back(energyExpression.differentiate(param).createCompiledExpression());
    }
    ixn = new ReferenceCustomCentroidBondIxn(force.getNumGroupsPerBond(), groupAtoms, normalizedWeights, bondGroups, energyExpression, bondParameterNames, distances, angles, dihedrals, energyParamDerivExpressions);
1752
1753
1754

    // Delete the custom functions.

peastman's avatar
peastman committed
1755
1756
    for (auto& function : functions)
        delete function.second;
1757
1758
1759
}

double ReferenceCalcCustomCentroidBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1760
1761
1762
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
1763
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1764
1765
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1766
1767
    if (usePeriodic)
        ixn->setPeriodic(extractBoxVectors(context));
1768
1769
1770
1771
1772
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
    ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
    return energy;
}

void ReferenceCalcCustomCentroidBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force) {
    if (numBonds != force.getNumBonds())
        throw OpenMMException("updateParametersInContext: The number of bonds has changed");

    // Record the values.

    int numParameters = force.getNumPerBondParameters();
    const vector<vector<int> >& bondGroups = ixn->getBondGroups();
    vector<int> groups;
    vector<double> params;
    for (int i = 0; i < numBonds; ++i) {
        force.getBondParameters(i, groups, params);
        for (int j = 0; j < groups.size(); j++)
            if (groups[j] != bondGroups[i][j])
                throw OpenMMException("updateParametersInContext: The set of groups in a bond has changed");
        for (int j = 0; j < numParameters; j++)
peastman's avatar
peastman committed
1792
            bondParamArray[i][j] = params[j];
1793
1794
1795
    }
}

1796
1797
1798
1799
1800
1801
ReferenceCalcCustomCompoundBondForceKernel::~ReferenceCalcCustomCompoundBondForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system, const CustomCompoundBondForce& force) {
1802
    usePeriodic = force.usesPeriodicBoundaryConditions();
1803
1804
1805
1806
1807
1808

    // Build the arrays.

    numBonds = force.getNumBonds();
    vector<vector<int> > bondParticles(numBonds);
    int numBondParameters = force.getNumPerBondParameters();
1809
1810
1811
    bondParamArray.resize(numBonds);
    for (int i = 0; i < numBonds; ++i)
        force.getBondParameters(i, bondParticles[i], bondParamArray[i]);
1812
1813
1814
1815

    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
1816
    for (int i = 0; i < force.getNumFunctions(); i++)
1817
        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829

    // Parse the expression and create the object used to calculate the interaction.

    map<string, vector<int> > distances;
    map<string, vector<int> > angles;
    map<string, vector<int> > dihedrals;
    Lepton::ParsedExpression energyExpression = CustomCompoundBondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
    vector<string> bondParameterNames;
    for (int i = 0; i < numBondParameters; i++)
        bondParameterNames.push_back(force.getPerBondParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
1830
1831
1832
1833
1834
1835
1836
    vector<Lepton::CompiledExpression> energyParamDerivExpressions;
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
        string param = force.getEnergyParameterDerivativeName(i);
        energyParamDerivNames.push_back(param);
        energyParamDerivExpressions.push_back(energyExpression.differentiate(param).createCompiledExpression());
    }
    ixn = new ReferenceCustomCompoundBondIxn(force.getNumParticlesPerBond(), bondParticles, energyExpression, bondParameterNames, distances, angles, dihedrals, energyParamDerivExpressions);
1837
1838
1839

    // Delete the custom functions.

peastman's avatar
peastman committed
1840
1841
    for (auto& function : functions)
        delete function.second;
1842
1843
1844
}

double ReferenceCalcCustomCompoundBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1845
1846
1847
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
1848
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1849
1850
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1851
1852
    if (usePeriodic)
        ixn->setPeriodic(extractBoxVectors(context));
1853
1854
1855
1856
1857
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
    ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
1858
1859
1860
    return energy;
}

1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
void ReferenceCalcCustomCompoundBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force) {
    if (numBonds != force.getNumBonds())
        throw OpenMMException("updateParametersInContext: The number of bonds has changed");

    // Record the values.

    int numParameters = force.getNumPerBondParameters();
    const vector<vector<int> >& bondAtoms = ixn->getBondAtoms();
    vector<int> particles;
    vector<double> params;
    for (int i = 0; i < numBonds; ++i) {
        force.getBondParameters(i, particles, params);
1873
        for (int j = 0; j < particles.size(); j++)
1874
1875
1876
            if (particles[j] != bondAtoms[i][j])
                throw OpenMMException("updateParametersInContext: The set of particles in a bond has changed");
        for (int j = 0; j < numParameters; j++)
peastman's avatar
peastman committed
1877
            bondParamArray[i][j] = params[j];
1878
1879
1880
    }
}

1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
ReferenceCalcCustomManyParticleForceKernel::~ReferenceCalcCustomManyParticleForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

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

    // Build the arrays.

    numParticles = system.getNumParticles();
    int numParticleParameters = force.getNumPerParticleParameters();
1892
    particleParamArray.resize(numParticles);
1893
1894
    for (int i = 0; i < numParticles; ++i) {
        int type;
1895
        force.getParticleParameters(i, particleParamArray[i], type);
1896
1897
1898
    }
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
1899
    ixn = new ReferenceCustomManyParticleIxn(force);
1900
1901
1902
1903
1904
    nonbondedMethod = CalcCustomManyParticleForceKernel::NonbondedMethod(force.getNonbondedMethod());
    cutoffDistance = force.getCutoffDistance();
}

double ReferenceCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1905
1906
1907
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
1908
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1909
1910
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1911
    if (nonbondedMethod == CutoffPeriodic) {
peastman's avatar
peastman committed
1912
        Vec3* boxVectors = extractBoxVectors(context);
1913
        double minAllowedSize = 2*cutoffDistance;
1914
        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
1915
            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
1916
        ixn->setPeriodic(boxVectors);
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
    }
    ixn->calculateIxn(posData, particleParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
    return energy;
}

void ReferenceCalcCustomManyParticleForceKernel::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
1935
            particleParamArray[i][j] = parameters[j];
1936
1937
1938
    }
}

1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
ReferenceCalcGayBerneForceKernel::~ReferenceCalcGayBerneForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

void ReferenceCalcGayBerneForceKernel::initialize(const System& system, const GayBerneForce& force) {
    ixn = new ReferenceGayBerneForce(force);
}

double ReferenceCalcGayBerneForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
    return ixn->calculateForce(extractPositions(context), extractForces(context), extractBoxVectors(context));
}

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

1958
1959
1960
1961
1962
ReferenceCalcCustomCVForceKernel::~ReferenceCalcCustomCVForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

1963
void ReferenceCalcCustomCVForceKernel::initialize(const System& system, const CustomCVForce& force, ContextImpl& innerContext) {
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++)
        energyParamDerivNames.push_back(force.getEnergyParameterDerivativeName(i));
    ixn = new ReferenceCustomCVForce(force);
}

double ReferenceCalcCustomCVForceKernel::execute(ContextImpl& context, ContextImpl& innerContext, bool includeForces, bool includeEnergy) {
    copyState(context, innerContext);
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
    map<string, double> globalParameters;
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    ixn->calculateIxn(innerContext, posData, globalParameters, forceData, includeEnergy ? &energy : NULL, energyParamDerivs);
    return energy;
}

void ReferenceCalcCustomCVForceKernel::copyState(ContextImpl& context, ContextImpl& innerContext) {
    extractPositions(innerContext) = extractPositions(context);
    extractVelocities(innerContext) = extractVelocities(context);
1987
1988
1989
1990
    Vec3 a, b, c;
    context.getPeriodicBoxVectors(a, b, c);
    innerContext.setPeriodicBoxVectors(a, b, c);
    innerContext.setTime(context.getTime());
1991
1992
1993
1994
1995
    map<string, double> innerParameters = innerContext.getParameters();
    for (auto& param : innerParameters)
        innerContext.setParameter(param.first, context.getParameter(param.first));
}

1996
1997
1998
1999
void ReferenceCalcCustomCVForceKernel::copyParametersToContext(ContextImpl& context, const CustomCVForce& force) {
    ixn->rebuildExpressions(force);
}

2000
2001
void ReferenceCalcRMSDForceKernel::initialize(const System& system, const RMSDForce& force) {
    particles = force.getParticles();
peastman's avatar
peastman committed
2002
2003
2004
    if (particles.size() == 0)
        for (int i = 0; i < system.getNumParticles(); i++)
            particles.push_back(i);
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
    referencePos = force.getReferencePositions();
    Vec3 center;
    for (int i : particles)
        center += referencePos[i];
    center /= particles.size();
    for (Vec3& p : referencePos)
        p -= center;
}

double ReferenceCalcRMSDForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    ReferenceRMSDForce rmsd(referencePos, particles);
    return rmsd.calculateIxn(posData, forceData);
}

void ReferenceCalcRMSDForceKernel::copyParametersToContext(ContextImpl& context, const RMSDForce& force) {
    if (referencePos.size() != force.getReferencePositions().size())
        throw OpenMMException("updateParametersInContext: The number of reference positions has changed");
    particles = force.getParticles();
peastman's avatar
peastman committed
2025
2026
2027
    if (particles.size() == 0)
        for (int i = 0; i < referencePos.size(); i++)
            particles.push_back(i);
2028
2029
2030
2031
2032
2033
2034
2035
2036
    referencePos = force.getReferencePositions();
    Vec3 center;
    for (int i : particles)
        center += referencePos[i];
    center /= particles.size();
    for (Vec3& p : referencePos)
        p -= center;
}

2037
2038
2039
2040
2041
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
    if (dynamics)
        delete dynamics;
}

2042
void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
2043
    int numParticles = system.getNumParticles();
2044
    masses.resize(numParticles);
2045
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2046
        masses[i] = system.getParticleMass(i);
2047
2048
}

2049
void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
2050
    double stepSize = integrator.getStepSize();
peastman's avatar
peastman committed
2051
2052
2053
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2054
2055
2056
    if (dynamics == 0 || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
2057
        if (dynamics)
2058
            delete dynamics;
peastman's avatar
peastman committed
2059
        dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), stepSize);
2060
        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
2061
2062
        prevStepSize = stepSize;
    }
2063
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
2064
    data.time += stepSize;
2065
    data.stepCount++;
2066
}
2067

2068
double ReferenceIntegrateVerletStepKernel::computeKineticEnergy(ContextImpl& context, const VerletIntegrator& integrator) {
2069
    return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
2070
2071
}

2072
2073
2074
2075
ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
}
2076

2077
void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
2078
    int numParticles = system.getNumParticles();
2079
    masses.resize(numParticles);
2080
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2081
        masses[i] = system.getParticleMass(i);
2082
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
2083
2084
}

2085
void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
2086
2087
2088
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
peastman's avatar
peastman committed
2089
2090
2091
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2092
2093
2094
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
2095
        if (dynamics)
2096
            delete dynamics;
2097
        dynamics = new ReferenceStochasticDynamics(
2098
                context.getSystem().getNumParticles(), 
peastman's avatar
peastman committed
2099
2100
2101
                stepSize, 
                friction, 
                temperature);
2102
        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
2103
2104
2105
2106
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
2107
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
2108
    data.time += stepSize;
2109
    data.stepCount++;
2110
2111
}

2112
double ReferenceIntegrateLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const LangevinIntegrator& integrator) {
2113
    return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
2114
2115
}

2116
2117
2118
2119
2120
ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
    if (dynamics)
        delete dynamics;
}

2121
void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
2122
    int numParticles = system.getNumParticles();
2123
    masses.resize(numParticles);
2124
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2125
        masses[i] = system.getParticleMass(i);
2126
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
2127
2128
}

2129
void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
2130
2131
2132
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
peastman's avatar
peastman committed
2133
2134
2135
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2136
2137
2138
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
2139
        if (dynamics)
2140
            delete dynamics;
2141
        dynamics = new ReferenceBrownianDynamics(
2142
                context.getSystem().getNumParticles(), 
peastman's avatar
peastman committed
2143
2144
2145
                stepSize, 
                friction, 
                temperature);
2146
        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
2147
2148
2149
2150
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
2151
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
2152
    data.time += stepSize;
2153
    data.stepCount++;
2154
2155
}

2156
double ReferenceIntegrateBrownianStepKernel::computeKineticEnergy(ContextImpl& context, const BrownianIntegrator& integrator) {
2157
    return computeShiftedKineticEnergy(context, masses, 0);
2158
2159
}

2160
2161
2162
2163
2164
2165
2166
ReferenceIntegrateVariableLangevinStepKernel::~ReferenceIntegrateVariableLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
}

void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& system, const VariableLangevinIntegrator& integrator) {
    int numParticles = system.getNumParticles();
2167
    masses.resize(numParticles);
2168
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2169
        masses[i] = system.getParticleMass(i);
2170
2171
2172
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
}

2173
double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
2174
2175
2176
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double errorTol = integrator.getErrorTolerance();
peastman's avatar
peastman committed
2177
2178
2179
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2180
2181
2182
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || errorTol != prevErrorTol) {
        // Recreate the computation objects with the new parameters.

2183
        if (dynamics)
2184
            delete dynamics;
peastman's avatar
peastman committed
2185
        dynamics = new ReferenceVariableStochasticDynamics(context.getSystem().getNumParticles(), friction, temperature, errorTol);
2186
        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
2187
2188
2189
2190
        prevTemp = temperature;
        prevFriction = friction;
        prevErrorTol = errorTol;
    }
peastman's avatar
peastman committed
2191
    double maxStepSize = maxTime-data.time;
2192
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize, integrator.getConstraintTolerance());
2193
2194
2195
2196
    data.time += dynamics->getDeltaT();
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
2197
    return dynamics->getDeltaT();
2198
2199
}

2200
double ReferenceIntegrateVariableLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const VariableLangevinIntegrator& integrator) {
2201
    return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
2202
2203
}

2204
2205
2206
2207
2208
2209
2210
ReferenceIntegrateVariableVerletStepKernel::~ReferenceIntegrateVariableVerletStepKernel() {
    if (dynamics)
        delete dynamics;
}

void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system, const VariableVerletIntegrator& integrator) {
    int numParticles = system.getNumParticles();
2211
    masses.resize(numParticles);
2212
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2213
        masses[i] = system.getParticleMass(i);
2214
2215
}

2216
double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
2217
    double errorTol = integrator.getErrorTolerance();
peastman's avatar
peastman committed
2218
2219
2220
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2221
    if (dynamics == 0 || errorTol != prevErrorTol) {
2222
2223
        // Recreate the computation objects with the new parameters.

2224
        if (dynamics)
2225
            delete dynamics;
peastman's avatar
peastman committed
2226
        dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), errorTol);
2227
        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
2228
        prevErrorTol = errorTol;
2229
    }
peastman's avatar
peastman committed
2230
    double maxStepSize = maxTime-data.time;
2231
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize, integrator.getConstraintTolerance());
2232
    data.time += dynamics->getDeltaT();
2233
2234
2235
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
2236
    return dynamics->getDeltaT();
2237
2238
}

2239
double ReferenceIntegrateVariableVerletStepKernel::computeKineticEnergy(ContextImpl& context, const VariableVerletIntegrator& integrator) {
2240
    return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
2241
2242
}

2243
2244
2245
2246
2247
2248
2249
2250
2251
ReferenceIntegrateCustomStepKernel::~ReferenceIntegrateCustomStepKernel() {
    if (dynamics)
        delete dynamics;
}

void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const CustomIntegrator& integrator) {
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2252
        masses[i] = system.getParticleMass(i);
2253
    perDofValues.resize(integrator.getNumPerDofVariables());
peastman's avatar
peastman committed
2254
2255
    for (auto& values : perDofValues)
        values.resize(numParticles);
2256
2257
2258
2259

    // Create the computation objects.

    dynamics = new ReferenceCustomDynamics(system.getNumParticles(), integrator);
2260
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
2261
2262
2263
}

void ReferenceIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid) {
peastman's avatar
peastman committed
2264
2265
2266
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2267
2268
2269
2270
2271
2272
2273
2274
2275
2276
    
    // Record global variables.
    
    map<string, double> globals;
    globals["dt"] = integrator.getStepSize();
    for (int i = 0; i < integrator.getNumGlobalVariables(); i++)
        globals[integrator.getGlobalVariableName(i)] = globalValues[i];
    
    // Execute the step.
    
2277
2278
    dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
    dynamics->update(context, context.getSystem().getNumParticles(), posData, velData, forceData, masses, globals, perDofValues, forcesAreValid, integrator.getConstraintTolerance());
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288
    
    // Record changed global variables.
    
    integrator.setStepSize(globals["dt"]);
    for (int i = 0; i < (int) globalValues.size(); i++)
        globalValues[i] = globals[integrator.getGlobalVariableName(i)];
    data.time += dynamics->getDeltaT();
    data.stepCount++;
}

2289
double ReferenceIntegrateCustomStepKernel::computeKineticEnergy(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid) {
peastman's avatar
peastman committed
2290
2291
2292
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2293
2294
2295
2296
2297
2298
2299
2300
2301
2302
2303
2304
2305
    
    // Record global variables.
    
    map<string, double> globals;
    globals["dt"] = integrator.getStepSize();
    for (int i = 0; i < integrator.getNumGlobalVariables(); i++)
        globals[integrator.getGlobalVariableName(i)] = globalValues[i];
    
    // Compute the kinetic energy.
    
    return dynamics->computeKineticEnergy(context, context.getSystem().getNumParticles(), posData, velData, forceData, masses, globals, perDofValues, forcesAreValid);
}

2306
2307
2308
2309
2310
2311
2312
2313
2314
2315
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
void ReferenceIntegrateCustomStepKernel::getGlobalVariables(ContextImpl& context, vector<double>& values) const {
    values = globalValues;
}

void ReferenceIntegrateCustomStepKernel::setGlobalVariables(ContextImpl& context, const vector<double>& values) {
    globalValues = values;
}

void ReferenceIntegrateCustomStepKernel::getPerDofVariable(ContextImpl& context, int variable, vector<Vec3>& values) const {
    values.resize(perDofValues[variable].size());
    for (int i = 0; i < (int) values.size(); i++)
        values[i] = perDofValues[variable][i];
}

void ReferenceIntegrateCustomStepKernel::setPerDofVariable(ContextImpl& context, int variable, const vector<Vec3>& values) {
    perDofValues[variable].resize(values.size());
    for (int i = 0; i < (int) values.size(); i++)
        perDofValues[variable][i] = values[i];
}

2326
2327
2328
2329
2330
ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
    if (thermostat)
        delete thermostat;
}

2331
void ReferenceApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
Peter Eastman's avatar
Peter Eastman committed
2332
    int numParticles = system.getNumParticles();
2333
    masses.resize(numParticles);
2334
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2335
        masses[i] = system.getParticleMass(i);
2336
    this->thermostat = new ReferenceAndersenThermostat();
2337
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
2338
    particleGroups = AndersenThermostatImpl::calcParticleGroups(system);
2339
2340
}

2341
void ReferenceApplyAndersenThermostatKernel::execute(ContextImpl& context) {
peastman's avatar
peastman committed
2342
    vector<Vec3>& velData = extractVelocities(context);
2343
    thermostat->applyThermostat(particleGroups, velData, masses,
peastman's avatar
peastman committed
2344
2345
2346
        context.getParameter(AndersenThermostat::Temperature()),
        context.getParameter(AndersenThermostat::CollisionFrequency()),
        context.getIntegrator().getStepSize());
2347
2348
}

2349
2350
2351
2352
2353
ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel() {
    if (barostat)
        delete barostat;
}

2354
void ReferenceApplyMonteCarloBarostatKernel::initialize(const System& system, const Force& barostat) {
2355
2356
}

2357
void ReferenceApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ) {
2358
2359
    if (barostat == NULL)
        barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), context.getMolecules());
peastman's avatar
peastman committed
2360
2361
    vector<Vec3>& posData = extractPositions(context);
    Vec3* boxVectors = extractBoxVectors(context);
2362
    barostat->applyBarostat(posData, boxVectors, scaleX, scaleY, scaleZ);
2363
2364
2365
}

void ReferenceApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
peastman's avatar
peastman committed
2366
    vector<Vec3>& posData = extractPositions(context);
2367
2368
2369
    barostat->restorePositions(posData);
}

2370
2371
void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
    frequency = force.getFrequency();
Peter Eastman's avatar
Peter Eastman committed
2372
    masses.resize(system.getNumParticles());
2373
    for (size_t i = 0; i < masses.size(); ++i)
Peter Eastman's avatar
Peter Eastman committed
2374
        masses[i] = system.getParticleMass(i);
2375
2376
}

2377
void ReferenceRemoveCMMotionKernel::execute(ContextImpl& context) {
2378
    if (data.stepCount%frequency != 0)
2379
        return;
peastman's avatar
peastman committed
2380
    vector<Vec3>& velData = extractVelocities(context);
2381
2382
2383
    
    // Calculate the center of mass momentum.
    
peastman's avatar
peastman committed
2384
2385
    double momentum[] = {0.0, 0.0, 0.0};
    double mass = 0.0;
2386
    for (size_t i = 0; i < masses.size(); ++i) {
peastman's avatar
peastman committed
2387
2388
2389
2390
        momentum[0] += masses[i]*velData[i][0];
        momentum[1] += masses[i]*velData[i][1];
        momentum[2] += masses[i]*velData[i][2];
        mass += masses[i];
2391
2392
    }
    
Peter Eastman's avatar
Peter Eastman committed
2393
    // Adjust the particle velocities.
2394
    
2395
2396
2397
    momentum[0] /= mass;
    momentum[1] /= mass;
    momentum[2] /= mass;
2398
    for (size_t i = 0; i < masses.size(); ++i) {
2399
2400
2401
2402
2403
        if (masses[i] != 0.0) {
            velData[i][0] -= momentum[0];
            velData[i][1] -= momentum[1];
            velData[i][2] -= momentum[2];
        }
2404
2405
    }
}