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
389
390
ReferenceCalcCustomBondForceKernel::~ReferenceCalcCustomBondForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

391
392
393
void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const CustomBondForce& force) {
    numBonds = force.getNumBonds();
    int numParameters = force.getNumPerBondParameters();
394
    usePeriodic = force.usesPeriodicBoundaryConditions();
395
396
397

    // Build the arrays.

398
399
    bondIndexArray.resize(numBonds, vector<int>(2));
    bondParamArray.resize(numBonds, vector<double>(numParameters));
400
    vector<double> params;
401
    for (int i = 0; i < numBonds; ++i) {
402
403
404
405
406
        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
407
            bondParamArray[i][j] = params[j];
408
409
410
411
412
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
413
414
    energyExpression = expression.createCompiledExpression();
    forceExpression = expression.differentiate("r").createCompiledExpression();
415
416
417
418
    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));
419
420
421
422
423
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
        string param = force.getEnergyParameterDerivativeName(i);
        energyParamDerivNames.push_back(param);
        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
    }
424
425
426
427
428
    set<string> variables;
    variables.insert("r");
    variables.insert(parameterNames.begin(), parameterNames.end());
    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
    validateVariables(expression.getRootNode(), variables);
429
    ixn = new ReferenceCustomBondIxn(energyExpression, forceExpression, parameterNames, energyParamDerivExpressions);
430
431
}

432
double ReferenceCalcCustomBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
433
434
435
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
436
    map<string, double> globalParameters;
peastman's avatar
peastman committed
437
438
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
439
    ixn->setGlobalParameters(globalParameters);
440
    if (usePeriodic)
441
        ixn->setPeriodic(extractBoxVectors(context));
442
443
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
    for (int i = 0; i < numBonds; i++)
444
        ixn->calculateBondIxn(bondIndexArray[i], posData, bondParamArray[i], forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
445
446
447
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
448
449
450
    return energy;
}

451
452
453
454
455
456
457
458
459
460
461
462
463
464
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
465
            bondParamArray[i][j] = params[j];
466
467
468
    }
}

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

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

498
499
500
501
502
503
504
505
506
507
508
509
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
510
511
        angleParamArray[i][0] = angle;
        angleParamArray[i][1] = k;
512
513
514
    }
}

515
516
517
518
519
ReferenceCalcCustomAngleForceKernel::~ReferenceCalcCustomAngleForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

520
521
522
void ReferenceCalcCustomAngleForceKernel::initialize(const System& system, const CustomAngleForce& force) {
    numAngles = force.getNumAngles();
    int numParameters = force.getNumPerAngleParameters();
523
    usePeriodic = force.usesPeriodicBoundaryConditions();
524
525
526

    // Build the arrays.

527
528
    angleIndexArray.resize(numAngles, vector<int>(3));
    angleParamArray.resize(numAngles, vector<double>(numParameters));
529
    vector<double> params;
530
    for (int i = 0; i < numAngles; ++i) {
531
532
533
534
535
536
        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
537
            angleParamArray[i][j] = params[j];
538
539
540
541
542
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
543
544
    energyExpression = expression.createCompiledExpression();
    forceExpression = expression.differentiate("theta").createCompiledExpression();
545
546
547
548
    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));
549
550
551
552
553
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
        string param = force.getEnergyParameterDerivativeName(i);
        energyParamDerivNames.push_back(param);
        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
    }
554
555
556
557
558
    set<string> variables;
    variables.insert("theta");
    variables.insert(parameterNames.begin(), parameterNames.end());
    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
    validateVariables(expression.getRootNode(), variables);
559
    ixn = new ReferenceCustomAngleIxn(energyExpression, forceExpression, parameterNames, energyParamDerivExpressions);
560
561
}

562
double ReferenceCalcCustomAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
563
564
565
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
566
    map<string, double> globalParameters;
peastman's avatar
peastman committed
567
568
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
569
    ixn->setGlobalParameters(globalParameters);
570
    if (usePeriodic)
571
        ixn->setPeriodic(extractBoxVectors(context));
572
573
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
    for (int i = 0; i < numAngles; i++)
574
        ixn->calculateBondIxn(angleIndexArray[i], posData, angleParamArray[i], forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
575
576
577
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
578
579
580
    return energy;
}

581
582
583
584
585
586
587
588
589
590
591
592
593
594
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
595
            angleParamArray[i][j] = params[j];
596
597
598
    }
}

599
600
void ReferenceCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
    numTorsions = force.getNumTorsions();
601
602
    torsionIndexArray.resize(numTorsions, vector<int>(4));
    torsionParamArray.resize(numTorsions, vector<double>(3));
603
    for (int i = 0; i < numTorsions; ++i) {
Peter Eastman's avatar
Peter Eastman committed
604
        int particle1, particle2, particle3, particle4, periodicity;
605
        double phase, k;
Peter Eastman's avatar
Peter Eastman committed
606
607
608
609
610
        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
611
612
613
        torsionParamArray[i][0] = k;
        torsionParamArray[i][1] = phase;
        torsionParamArray[i][2] = periodicity;
614
    }
615
    usePeriodic = force.usesPeriodicBoundaryConditions();
616
617
}

618
double ReferenceCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
619
620
621
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
622
623
    ReferenceBondForce refBondForce;
    ReferenceProperDihedralBond periodicTorsionBond;
624
625
    if (usePeriodic)
        periodicTorsionBond.setPeriodic(extractBoxVectors(context));
626
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, periodicTorsionBond);
627
628
629
    return energy;
}

630
631
632
633
634
635
636
637
638
639
640
641
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
642
643
644
        torsionParamArray[i][0] = k;
        torsionParamArray[i][1] = phase;
        torsionParamArray[i][2] = periodicity;
645
646
647
    }
}

648
649
void ReferenceCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
    numTorsions = force.getNumTorsions();
650
651
    torsionIndexArray.resize(numTorsions, vector<int>(4));
    torsionParamArray.resize(numTorsions, vector<double>(6));
652
    for (int i = 0; i < numTorsions; ++i) {
Peter Eastman's avatar
Peter Eastman committed
653
        int particle1, particle2, particle3, particle4;
654
        double c0, c1, c2, c3, c4, c5;
Peter Eastman's avatar
Peter Eastman committed
655
656
657
658
659
        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
660
661
662
663
664
665
        torsionParamArray[i][0] = c0;
        torsionParamArray[i][1] = c1;
        torsionParamArray[i][2] = c2;
        torsionParamArray[i][3] = c3;
        torsionParamArray[i][4] = c4;
        torsionParamArray[i][5] = c5;
666
    }
667
    usePeriodic = force.usesPeriodicBoundaryConditions();
668
669
}

670
double ReferenceCalcRBTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
671
672
673
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
674
675
    ReferenceBondForce refBondForce;
    ReferenceRbDihedralBond rbTorsionBond;
676
677
    if (usePeriodic)
        rbTorsionBond.setPeriodic(extractBoxVectors(context));
678
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, rbTorsionBond);
679
680
681
    return energy;
}

682
683
684
685
686
687
688
689
690
691
692
693
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
694
695
696
697
698
699
        torsionParamArray[i][0] = c0;
        torsionParamArray[i][1] = c1;
        torsionParamArray[i][2] = c2;
        torsionParamArray[i][3] = c3;
        torsionParamArray[i][4] = c4;
        torsionParamArray[i][5] = c5;
700
701
702
    }
}

703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
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]);
    }
727
    usePeriodic = force.usesPeriodicBoundaryConditions();
728
729
}

730
double ReferenceCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
731
732
733
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double totalEnergy = 0;
734
    ReferenceCMAPTorsionIxn torsion(coeff, torsionMaps, torsionIndices);
735
736
    if (usePeriodic)
        torsion.setPeriodic(extractBoxVectors(context));
737
738
739
740
    torsion.calculateIxn(posData, forceData, &totalEnergy);
    return totalEnergy;
}

741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
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");
    }
}

775
776
777
778
779
ReferenceCalcCustomTorsionForceKernel::~ReferenceCalcCustomTorsionForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

780
781
782
void ReferenceCalcCustomTorsionForceKernel::initialize(const System& system, const CustomTorsionForce& force) {
    numTorsions = force.getNumTorsions();
    int numParameters = force.getNumPerTorsionParameters();
783
    usePeriodic = force.usesPeriodicBoundaryConditions();
784
785
786

    // Build the arrays.

787
788
    torsionIndexArray.resize(numTorsions, vector<int>(4));
    torsionParamArray.resize(numTorsions, vector<double>(numParameters));
789
    vector<double> params;
790
    for (int i = 0; i < numTorsions; ++i) {
791
792
793
794
795
796
797
        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
798
            torsionParamArray[i][j] = params[j];
799
800
801
802
803
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
804
805
    energyExpression = expression.createCompiledExpression();
    forceExpression = expression.differentiate("theta").createCompiledExpression();
806
807
808
809
    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));
810
811
812
813
814
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
        string param = force.getEnergyParameterDerivativeName(i);
        energyParamDerivNames.push_back(param);
        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
    }
815
816
817
818
819
    set<string> variables;
    variables.insert("theta");
    variables.insert(parameterNames.begin(), parameterNames.end());
    variables.insert(globalParameterNames.begin(), globalParameterNames.end());
    validateVariables(expression.getRootNode(), variables);
820
    ixn = new ReferenceCustomTorsionIxn(energyExpression, forceExpression, parameterNames, energyParamDerivExpressions);
821
822
}

823
double ReferenceCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
824
825
826
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
827
    map<string, double> globalParameters;
peastman's avatar
peastman committed
828
829
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
830
    ixn->setGlobalParameters(globalParameters);
831
    if (usePeriodic)
832
        ixn->setPeriodic(extractBoxVectors(context));
833
834
    vector<double> energyParamDerivValues(energyParamDerivNames.size()+1, 0.0);
    for (int i = 0; i < numTorsions; i++)
835
        ixn->calculateBondIxn(torsionIndexArray[i], posData, torsionParamArray[i], forceData, includeEnergy ? &energy : NULL, &energyParamDerivValues[0]);
836
837
838
    map<string, double>& energyParamDerivs = extractEnergyParameterDerivatives(context);
    for (int i = 0; i < energyParamDerivNames.size(); i++)
        energyParamDerivs[energyParamDerivNames[i]] += energyParamDerivValues[i];
839
840
841
    return energy;
}

842
843
844
845
846
847
848
849
850
851
852
853
854
855
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
856
            torsionParamArray[i][j] = params[j];
857
858
859
    }
}

860
861
862
863
864
ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
    if (neighborList != NULL)
        delete neighborList;
}

865
866
867
868
void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {

    // Identify which exceptions are 1-4 interactions.

Peter Eastman's avatar
Peter Eastman committed
869
    numParticles = force.getNumParticles();
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
    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();
885
886
887
    bonded14IndexArray.resize(num14, vector<int>(2));
    bonded14ParamArray.resize(num14, vector<double>(3));
    particleParamArray.resize(numParticles, vector<double>(3));
888
889
890
891
    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]);
892
    this->exclusions = exclusions;
893
    for (int i = 0; i < num14; ++i) {
Peter Eastman's avatar
Peter Eastman committed
894
        int particle1, particle2;
895
        force.getExceptionParameters(nb14s[i], particle1, particle2, baseExceptionParams[i][0], baseExceptionParams[i][1], baseExceptionParams[i][2]);
Peter Eastman's avatar
Peter Eastman committed
896
897
        bonded14IndexArray[i][0] = particle1;
        bonded14IndexArray[i][1] = particle2;
898
899
900
901
902
903
904
905
906
907
908
909
910
911
    }
    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};
912
    }
913
    nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
peastman's avatar
peastman committed
914
    nonbondedCutoff = force.getCutoffDistance();
915
    if (nonbondedMethod == NoCutoff) {
916
        neighborList = NULL;
917
918
919
        useSwitchingFunction = false;
    }
    else {
920
        neighborList = new NeighborList();
921
922
923
        useSwitchingFunction = force.getUseSwitchingFunction();
        switchingDistance = force.getSwitchingDistance();
    }
924
925
926
    if (nonbondedMethod == Ewald) {
        double alpha;
        NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmax[0], kmax[1], kmax[2]);
peastman's avatar
peastman committed
927
        ewaldAlpha = alpha;
928
929
930
    }
    else if (nonbondedMethod == PME) {
        double alpha;
931
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2], false);
peastman's avatar
peastman committed
932
        ewaldAlpha = alpha;
933
    }
934
935
936
    else if (nonbondedMethod == LJPME) {
        double alpha;
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2], false);
peastman's avatar
peastman committed
937
        ewaldAlpha = alpha;
938
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, dispersionGridSize[0], dispersionGridSize[1], dispersionGridSize[2], true);
peastman's avatar
peastman committed
939
        ewaldDispersionAlpha = alpha;
940
941
        useSwitchingFunction = false;
    }
peastman's avatar
peastman committed
942
    rfDielectric = force.getReactionFieldDielectric();
943
944
945
946
    if (force.getUseDispersionCorrection())
        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
    else
        dispersionCoefficient = 0.0;
947
948
}

949
double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
950
    computeParameters(context);
peastman's avatar
peastman committed
951
952
953
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
954
    ReferenceLJCoulombIxn clj;
955
    bool periodic = (nonbondedMethod == CutoffPeriodic);
956
    bool ewald  = (nonbondedMethod == Ewald);
957
    bool pme  = (nonbondedMethod == PME);
958
    bool ljpme = (nonbondedMethod == LJPME);
959
    if (nonbondedMethod != NoCutoff) {
960
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic || ewald || pme || ljpme, nonbondedCutoff, 0.0);
961
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
962
    }
963
    if (periodic || ewald || pme || ljpme) {
peastman's avatar
peastman committed
964
        Vec3* boxVectors = extractBoxVectors(context);
965
        double minAllowedSize = 1.999999*nonbondedCutoff;
966
        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
967
            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
968
        clj.setPeriodic(boxVectors);
969
    }
970
971
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
972
    if (pme)
973
        clj.setUsePME(ewaldAlpha, gridSize);
974
975
976
977
    if (ljpme){
        clj.setUsePME(ewaldAlpha, gridSize);
        clj.setUseLJPME(ewaldDispersionAlpha, dispersionGridSize);
    }
978
979
    if (useSwitchingFunction)
        clj.setUseSwitchingFunction(switchingDistance);
980
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusions, forceData, includeEnergy ? &energy : NULL, includeDirect, includeReciprocal);
981
982
983
984
    if (includeDirect) {
        ReferenceBondForce refBondForce;
        ReferenceLJCoulomb14 nonbonded14;
        refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
985
        if (periodic || ewald || pme) {
peastman's avatar
peastman committed
986
            Vec3* boxVectors = extractBoxVectors(context);
987
            energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
988
        }
989
    }
990
991
992
    return energy;
}

993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
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.

1009
1010
    for (int i = 0; i < numParticles; ++i)
        force.getParticleParameters(i, baseParticleParams[i][0], baseParticleParams[i][1], baseParticleParams[i][2]);
1011
1012
    for (int i = 0; i < num14; ++i) {
        int particle1, particle2;
1013
        force.getExceptionParameters(nb14s[i], particle1, particle2, baseExceptionParams[i][0], baseExceptionParams[i][1], baseExceptionParams[i][2]);
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
        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);
}

1025
void ReferenceCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
1026
1027
    if (nonbondedMethod != PME && nonbondedMethod != LJPME)
        throw OpenMMException("getPMEParametersInContext: This Context is not using PME or LJPME");
1028
1029
1030
1031
1032
1033
    alpha = ewaldAlpha;
    nx = gridSize[0];
    ny = gridSize[1];
    nz = gridSize[2];
}

1034
void ReferenceCalcNonbondedForceKernel::getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
1035
1036
    if (nonbondedMethod != LJPME)
        throw OpenMMException("getPMEParametersInContext: This Context is not using LJPME");
1037
1038
1039
1040
    alpha = ewaldDispersionAlpha;
    nx = dispersionGridSize[0];
    ny = dispersionGridSize[1];
    nz = dispersionGridSize[2];
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
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
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];
    }
}

1089
1090
1091
ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
    if (neighborList != NULL)
        delete neighborList;
1092
1093
    if (forceCopy != NULL)
        delete forceCopy;
1094
1095
1096
1097
}

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

1098
    // Record the exclusions.
1099
1100
1101

    numParticles = force.getNumParticles();
    exclusions.resize(numParticles);
1102
    for (int i = 0; i < force.getNumExclusions(); i++) {
1103
        int particle1, particle2;
1104
        force.getExclusionParticles(i, particle1, particle2);
1105
1106
1107
1108
1109
1110
        exclusions[particle1].insert(particle2);
        exclusions[particle2].insert(particle1);
    }

    // Build the arrays.

1111
    int numParameters = force.getNumPerParticleParameters();
1112
1113
1114
    particleParamArray.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        force.getParticleParameters(i, particleParamArray[i]);
1115
    nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
peastman's avatar
peastman committed
1116
    nonbondedCutoff = force.getCutoffDistance();
1117
    if (nonbondedMethod == NoCutoff) {
1118
        neighborList = NULL;
1119
1120
1121
        useSwitchingFunction = false;
    }
    else {
1122
        neighborList = new NeighborList();
1123
1124
1125
        useSwitchingFunction = force.getUseSwitchingFunction();
        switchingDistance = force.getSwitchingDistance();
    }
1126

1127
1128
1129
    // Create custom functions for the tabulated functions.

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

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

1135
    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
1136
1137
    energyExpression = expression.createCompiledExpression();
    forceExpression = expression.differentiate("r").createCompiledExpression();
1138
1139
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerParticleParameterName(i));
1140
    for (int i = 0; i < force.getNumGlobalParameters(); i++) {
1141
        globalParameterNames.push_back(force.getGlobalParameterName(i));
1142
1143
        globalParamValues[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
    }
1144
1145
1146
1147
1148
    for (int i = 0; i < force.getNumEnergyParameterDerivatives(); i++) {
        string param = force.getEnergyParameterDerivativeName(i);
        energyParamDerivNames.push_back(param);
        energyParamDerivExpressions.push_back(expression.differentiate(param).createCompiledExpression());
    }
1149
1150
1151
1152
1153
1154
1155
1156
    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);
1157
1158
1159

    // Delete the custom functions.

peastman's avatar
peastman committed
1160
1161
    for (auto& function : functions)
        delete function.second;
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
    
    // 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;
    }
1173
1174
1175
1176
1177
1178
1179
1180
    
    // 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));
    }
1181
1182
}

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

1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
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
1242
            particleParamArray[i][j] = parameters[j];
1243
    }
1244
1245
1246
1247
    
    // If necessary, recompute the long range correction.
    
    if (forceCopy != NULL) {
1248
        CustomNonbondedForceImpl::calcLongRangeCorrection(force, context.getOwner(), longRangeCoefficient, longRangeCoefficientDerivs);
1249
1250
1251
        hasInitializedLongRangeCorrection = true;
        *forceCopy = force;
    }
1252
1253
}

1254
ReferenceCalcGBSAOBCForceKernel::~ReferenceCalcGBSAOBCForceKernel() {
1255
    if (obc) {
Peter Eastman's avatar
Peter Eastman committed
1256
        delete obc->getObcParameters();
1257
1258
1259
1260
        delete obc;
    }
}

1261
void ReferenceCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
Peter Eastman's avatar
Peter Eastman committed
1262
1263
    int numParticles = system.getNumParticles();
    charges.resize(numParticles);
peastman's avatar
peastman committed
1264
1265
    vector<double> atomicRadii(numParticles);
    vector<double> scaleFactors(numParticles);
Peter Eastman's avatar
Peter Eastman committed
1266
    for (int i = 0; i < numParticles; ++i) {
1267
        double charge, radius, scalingFactor;
Peter Eastman's avatar
Peter Eastman committed
1268
        force.getParticleParameters(i, charge, radius, scalingFactor);
peastman's avatar
peastman committed
1269
1270
1271
        charges[i] = charge;
        atomicRadii[i] = radius;
        scaleFactors[i] = scalingFactor;
1272
    }
1273
    ObcParameters* obcParameters = new ObcParameters(numParticles, ObcParameters::ObcTypeII);
1274
    obcParameters->setAtomicRadii(atomicRadii);
1275
    obcParameters->setScaledRadiusFactors(scaleFactors);
peastman's avatar
peastman committed
1276
1277
    obcParameters->setSolventDielectric(force.getSolventDielectric());
    obcParameters->setSoluteDielectric(force.getSoluteDielectric());
1278
    obcParameters->setPi4Asolv(4*M_PI*force.getSurfaceAreaEnergy());
1279
    if (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff)
peastman's avatar
peastman committed
1280
        obcParameters->setUseCutoff(force.getCutoffDistance());
1281
    isPeriodic = (force.getNonbondedMethod() == GBSAOBCForce::CutoffPeriodic);
1282
    obc = new ReferenceObc(obcParameters);
1283
    obc->setIncludeAceApproximation(true);
1284
1285
}

1286
double ReferenceCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1287
1288
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
1289
    if (isPeriodic)
1290
        obc->getObcParameters()->setPeriodic(extractBoxVectors(context));
Mark Friedrichs's avatar
Mark Friedrichs committed
1291
    return obc->computeBornEnergyForces(posData, charges, forceData);
1292
1293
}

1294
1295
1296
1297
1298
1299
1300
1301
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
1302
1303
    vector<double> atomicRadii(numParticles);
    vector<double> scaleFactors(numParticles);
1304
1305
1306
    for (int i = 0; i < numParticles; ++i) {
        double charge, radius, scalingFactor;
        force.getParticleParameters(i, charge, radius, scalingFactor);
peastman's avatar
peastman committed
1307
1308
1309
        charges[i] = charge;
        atomicRadii[i] = radius;
        scaleFactors[i] = scalingFactor;
1310
1311
1312
1313
1314
    }
    obcParameters->setAtomicRadii(atomicRadii);
    obcParameters->setScaledRadiusFactors(scaleFactors);
}

1315
1316
1317
1318
1319
1320
ReferenceCalcCustomGBForceKernel::~ReferenceCalcCustomGBForceKernel() {
    if (neighborList != NULL)
        delete neighborList;
}

void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
    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.");
        }
    }
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347

    // 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();
1348
1349
1350
    particleParamArray.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        force.getParticleParameters(i, particleParamArray[i]);
1351
1352
1353
1354
1355
    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
1356
    nonbondedCutoff = force.getCutoffDistance();
1357
1358
1359
1360
1361
1362
1363
1364
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();

    // Create custom functions for the tabulated functions.

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

    // Parse the expressions for computed values.

1370
    valueDerivExpressions.resize(force.getNumComputedValues());
1371
    valueGradientExpressions.resize(force.getNumComputedValues());
1372
    valueParamDerivExpressions.resize(force.getNumComputedValues());
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
    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());
1385
1386
1387
1388
1389
    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();
1390
        valueExpressions.push_back(ex.createCompiledExpression());
1391
1392
        valueTypes.push_back(type);
        valueNames.push_back(name);
1393
        if (i == 0) {
1394
            valueDerivExpressions[i].push_back(ex.differentiate("r").createCompiledExpression());
1395
1396
            validateVariables(ex.getRootNode(), pairVariables);
        }
1397
        else {
1398
1399
1400
            valueGradientExpressions[i].push_back(ex.differentiate("x").createCompiledExpression());
            valueGradientExpressions[i].push_back(ex.differentiate("y").createCompiledExpression());
            valueGradientExpressions[i].push_back(ex.differentiate("z").createCompiledExpression());
1401
            for (int j = 0; j < i; j++)
1402
                valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).createCompiledExpression());
1403
            validateVariables(ex.getRootNode(), particleVariables);
1404
        }
1405
1406
1407
1408
1409
        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());
        }
1410
1411
1412
        particleVariables.insert(name);
        pairVariables.insert(name+"1");
        pairVariables.insert(name+"2");
1413
1414
    }

1415
    // Parse the expressions for energy terms.
1416
1417

    energyDerivExpressions.resize(force.getNumEnergyTerms());
1418
    energyGradientExpressions.resize(force.getNumEnergyTerms());
1419
    energyParamDerivExpressions.resize(force.getNumEnergyTerms());
1420
1421
1422
1423
1424
    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();
1425
        energyExpressions.push_back(ex.createCompiledExpression());
1426
1427
        energyTypes.push_back(type);
        if (type != CustomGBForce::SingleParticle)
1428
            energyDerivExpressions[i].push_back(ex.differentiate("r").createCompiledExpression());
1429
        for (int j = 0; j < force.getNumComputedValues(); j++) {
1430
            if (type == CustomGBForce::SingleParticle) {
1431
1432
1433
1434
                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());
1435
                validateVariables(ex.getRootNode(), particleVariables);
1436
            }
1437
            else {
1438
1439
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").createCompiledExpression());
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").createCompiledExpression());
1440
                validateVariables(ex.getRootNode(), pairVariables);
1441
1442
            }
        }
1443
1444
        for (int j = 0; j < force.getNumEnergyParameterDerivatives(); j++)
            energyParamDerivExpressions[i].push_back(ex.differentiate(force.getEnergyParameterDerivativeName(j)).createCompiledExpression());
1445
1446
1447
1448
    }

    // Delete the custom functions.

peastman's avatar
peastman committed
1449
1450
    for (auto& function : functions)
        delete function.second;
1451
1452
}

1453
double ReferenceCalcCustomGBForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1454
1455
1456
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
1457
1458
    ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueParamDerivExpressions, valueNames, valueTypes,
        energyExpressions, energyDerivExpressions, energyGradientExpressions, energyParamDerivExpressions, energyTypes, particleParameterNames);
1459
    bool periodic = (nonbondedMethod == CutoffPeriodic);
1460
    if (periodic)
1461
        ixn.setPeriodic(extractBoxVectors(context));
1462
    if (nonbondedMethod != NoCutoff) {
Peter Eastman's avatar
Peter Eastman committed
1463
1464
        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);
1465
1466
1467
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1468
1469
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1470
1471
1472
1473
1474
    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];
1475
1476
1477
    return energy;
}

1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
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
1490
            particleParamArray[i][j] = parameters[j];
1491
1492
1493
    }
}

peastman's avatar
peastman committed
1494
ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::PeriodicDistanceFunction(Vec3** boxVectorHandle) : boxVectorHandle(boxVectorHandle) {
1495
1496
1497
1498
1499
1500
1501
}

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

double ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::evaluate(const double* arguments) const {
peastman's avatar
peastman committed
1502
1503
    Vec3* boxVectors = *boxVectorHandle;
    Vec3 delta = Vec3(arguments[0], arguments[1], arguments[2])-Vec3(arguments[3], arguments[4], arguments[5]);
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
    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
1519
1520
    Vec3* boxVectors = *boxVectorHandle;
    Vec3 delta = Vec3(arguments[0], arguments[1], arguments[2])-Vec3(arguments[3], arguments[4], arguments[5]);
1521
1522
1523
1524
    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));
1525
1526
    if (r == 0)
        return 0.0;    
1527
1528
1529
1530
1531
1532
1533
1534
1535
    if (argIndex < 3)
        return delta[argIndex]/r;
    return -delta[argIndex-3]/r;
}

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

1536
1537
1538
1539
1540
ReferenceCalcCustomExternalForceKernel::~ReferenceCalcCustomExternalForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

1541
1542
1543
1544
1545
1546
1547
void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, const CustomExternalForce& force) {
    numParticles = force.getNumParticles();
    int numParameters = force.getNumPerParticleParameters();

    // Build the arrays.

    particles.resize(numParticles);
1548
1549
1550
    particleParamArray.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        force.getParticleParameters(i, particles[i], particleParamArray[i]);
1551
1552
1553

    // Parse the expression used to calculate the force.

1554
1555
1556
1557
    map<string, Lepton::CustomFunction*> functions;
    PeriodicDistanceFunction periodicDistance(&boxVectors);
    functions["periodicdistance"] = &periodicDistance;
    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
1558
1559
1560
1561
    energyExpression = expression.createCompiledExpression();
    forceExpressionX = expression.differentiate("x").createCompiledExpression();
    forceExpressionY = expression.differentiate("y").createCompiledExpression();
    forceExpressionZ = expression.differentiate("z").createCompiledExpression();
1562
1563
1564
1565
    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));
1566
1567
1568
1569
1570
1571
1572
    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);
1573
1574
    ixn = new ReferenceCustomExternalIxn(energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ, parameterNames);

1575
1576
}

1577
double ReferenceCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1578
1579
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
1580
    boxVectors = extractBoxVectors(context);
peastman's avatar
peastman committed
1581
    double energy = 0;
1582
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1583
1584
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1585
    ixn->setGlobalParameters(globalParameters);
1586
    for (int i = 0; i < numParticles; ++i)
1587
        ixn->calculateForce(particles[i], posData, particleParamArray[i], forceData, includeEnergy ? &energy : NULL);
1588
1589
1590
    return energy;
}

1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
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
1606
            particleParamArray[i][j] = parameters[j];
1607
1608
1609
    }
}

1610
ReferenceCalcCustomHbondForceKernel::~ReferenceCalcCustomHbondForceKernel() {
1611
1612
    if (ixn != NULL)
        delete ixn;
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
}

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.

1631
    vector<vector<int> > donorParticles(numDonors);
1632
    int numDonorParameters = force.getNumPerDonorParameters();
1633
    donorParamArray.resize(numDonors);
1634
    for (int i = 0; i < numDonors; ++i) {
1635
        int d1, d2, d3;
1636
        force.getDonorParameters(i, d1, d2, d3, donorParamArray[i]);
1637
1638
1639
        donorParticles[i].push_back(d1);
        donorParticles[i].push_back(d2);
        donorParticles[i].push_back(d3);
1640
    }
1641
    vector<vector<int> > acceptorParticles(numAcceptors);
1642
    int numAcceptorParameters = force.getNumPerAcceptorParameters();
1643
    acceptorParamArray.resize(numAcceptors);
1644
    for (int i = 0; i < numAcceptors; ++i) {
1645
        int a1, a2, a3;
1646
        force.getAcceptorParameters(i, a1, a2, a3, acceptorParamArray[i]);
1647
1648
1649
        acceptorParticles[i].push_back(a1);
        acceptorParticles[i].push_back(a2);
        acceptorParticles[i].push_back(a3);
1650
    }
1651
    NonbondedMethod nonbondedMethod = CalcCustomHbondForceKernel::NonbondedMethod(force.getNonbondedMethod());
peastman's avatar
peastman committed
1652
    nonbondedCutoff = force.getCutoffDistance();
1653
1654
1655
1656

    // Create custom functions for the tabulated functions.

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

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

1662
1663
1664
    map<string, vector<int> > distances;
    map<string, vector<int> > angles;
    map<string, vector<int> > dihedrals;
1665
    Lepton::ParsedExpression energyExpression = CustomHbondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
1666
1667
    vector<string> donorParameterNames;
    vector<string> acceptorParameterNames;
1668
1669
1670
1671
1672
1673
    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));
1674
    ixn = new ReferenceCustomHbondIxn(donorParticles, acceptorParticles, energyExpression, donorParameterNames, acceptorParameterNames, distances, angles, dihedrals);
1675
    isPeriodic = (nonbondedMethod == CutoffPeriodic);
1676
1677
    if (nonbondedMethod != NoCutoff)
        ixn->setUseCutoff(nonbondedCutoff);
1678
1679
1680

    // Delete the custom functions.

peastman's avatar
peastman committed
1681
1682
    for (auto& function : functions)
        delete function.second;
1683
1684
}

1685
double ReferenceCalcCustomHbondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1686
1687
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
1688
    if (isPeriodic)
1689
        ixn->setPeriodic(extractBoxVectors(context));
peastman's avatar
peastman committed
1690
    double energy = 0;
1691
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1692
1693
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1694
    ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusions, globalParameters, forceData, includeEnergy ? &energy : NULL);
1695
1696
1697
    return energy;
}

1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
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
1715
            donorParamArray[i][j] = parameters[j];
1716
1717
1718
1719
1720
1721
1722
1723
1724
    }
    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
1725
            acceptorParamArray[i][j] = parameters[j];
1726
1727
1728
    }
}

1729
1730
1731
1732
1733
1734
ReferenceCalcCustomCentroidBondForceKernel::~ReferenceCalcCustomCentroidBondForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

void ReferenceCalcCustomCentroidBondForceKernel::initialize(const System& system, const CustomCentroidBondForce& force) {
1735
    usePeriodic = force.usesPeriodicBoundaryConditions();
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748

    // 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();
1749
1750
1751
    bondParamArray.resize(numBonds);
    for (int i = 0; i < numBonds; ++i)
        force.getBondParameters(i, bondGroups[i], bondParamArray[i]);
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769

    // 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));
1770
1771
1772
1773
1774
1775
1776
    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);
1777
1778
1779

    // Delete the custom functions.

peastman's avatar
peastman committed
1780
1781
    for (auto& function : functions)
        delete function.second;
1782
1783
1784
}

double ReferenceCalcCustomCentroidBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1785
1786
1787
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
1788
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1789
1790
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1791
1792
    if (usePeriodic)
        ixn->setPeriodic(extractBoxVectors(context));
1793
1794
1795
1796
1797
    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];
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
    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
1817
            bondParamArray[i][j] = params[j];
1818
1819
1820
    }
}

1821
1822
1823
1824
1825
1826
ReferenceCalcCustomCompoundBondForceKernel::~ReferenceCalcCustomCompoundBondForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system, const CustomCompoundBondForce& force) {
1827
    usePeriodic = force.usesPeriodicBoundaryConditions();
1828
1829
1830
1831
1832
1833

    // Build the arrays.

    numBonds = force.getNumBonds();
    vector<vector<int> > bondParticles(numBonds);
    int numBondParameters = force.getNumPerBondParameters();
1834
1835
1836
    bondParamArray.resize(numBonds);
    for (int i = 0; i < numBonds; ++i)
        force.getBondParameters(i, bondParticles[i], bondParamArray[i]);
1837
1838
1839
1840

    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
1841
    for (int i = 0; i < force.getNumFunctions(); i++)
1842
        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854

    // 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));
1855
1856
1857
1858
1859
1860
1861
    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);
1862
1863
1864

    // Delete the custom functions.

peastman's avatar
peastman committed
1865
1866
    for (auto& function : functions)
        delete function.second;
1867
1868
1869
}

double ReferenceCalcCustomCompoundBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1870
1871
1872
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
1873
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1874
1875
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1876
1877
    if (usePeriodic)
        ixn->setPeriodic(extractBoxVectors(context));
1878
1879
1880
1881
1882
    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];
1883
1884
1885
    return energy;
}

1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
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);
1898
        for (int j = 0; j < particles.size(); j++)
1899
1900
1901
            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
1902
            bondParamArray[i][j] = params[j];
1903
1904
1905
    }
}

1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
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();
1917
    particleParamArray.resize(numParticles);
1918
1919
    for (int i = 0; i < numParticles; ++i) {
        int type;
1920
        force.getParticleParameters(i, particleParamArray[i], type);
1921
1922
1923
    }
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
1924
    ixn = new ReferenceCustomManyParticleIxn(force);
1925
1926
1927
1928
1929
    nonbondedMethod = CalcCustomManyParticleForceKernel::NonbondedMethod(force.getNonbondedMethod());
    cutoffDistance = force.getCutoffDistance();
}

double ReferenceCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
1930
1931
1932
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = 0;
1933
    map<string, double> globalParameters;
peastman's avatar
peastman committed
1934
1935
    for (auto& name : globalParameterNames)
        globalParameters[name] = context.getParameter(name);
1936
    if (nonbondedMethod == CutoffPeriodic) {
peastman's avatar
peastman committed
1937
        Vec3* boxVectors = extractBoxVectors(context);
1938
        double minAllowedSize = 2*cutoffDistance;
1939
        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
1940
            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
1941
        ixn->setPeriodic(boxVectors);
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
    }
    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
1960
            particleParamArray[i][j] = parameters[j];
1961
1962
1963
    }
}

1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
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);
}

1983
1984
1985
1986
1987
ReferenceCalcCustomCVForceKernel::~ReferenceCalcCustomCVForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

1988
void ReferenceCalcCustomCVForceKernel::initialize(const System& system, const CustomCVForce& force, ContextImpl& innerContext) {
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
    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);
2012
2013
2014
2015
    Vec3 a, b, c;
    context.getPeriodicBoxVectors(a, b, c);
    innerContext.setPeriodicBoxVectors(a, b, c);
    innerContext.setTime(context.getTime());
2016
2017
2018
2019
2020
    map<string, double> innerParameters = innerContext.getParameters();
    for (auto& param : innerParameters)
        innerContext.setParameter(param.first, context.getParameter(param.first));
}

2021
void ReferenceCalcCustomCVForceKernel::copyParametersToContext(ContextImpl& context, const CustomCVForce& force) {
2022
    ixn->updateTabulatedFunctions(force);
2023
2024
}

2025
2026
void ReferenceCalcRMSDForceKernel::initialize(const System& system, const RMSDForce& force) {
    particles = force.getParticles();
peastman's avatar
peastman committed
2027
2028
2029
    if (particles.size() == 0)
        for (int i = 0; i < system.getNumParticles(); i++)
            particles.push_back(i);
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
    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
2050
2051
2052
    if (particles.size() == 0)
        for (int i = 0; i < referencePos.size(); i++)
            particles.push_back(i);
2053
2054
2055
2056
2057
2058
2059
2060
2061
    referencePos = force.getReferencePositions();
    Vec3 center;
    for (int i : particles)
        center += referencePos[i];
    center /= particles.size();
    for (Vec3& p : referencePos)
        p -= center;
}

2062
2063
2064
2065
2066
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
    if (dynamics)
        delete dynamics;
}

2067
void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
2068
    int numParticles = system.getNumParticles();
2069
    masses.resize(numParticles);
2070
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2071
        masses[i] = system.getParticleMass(i);
2072
2073
}

2074
void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
2075
    double stepSize = integrator.getStepSize();
peastman's avatar
peastman committed
2076
2077
2078
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2079
2080
2081
    if (dynamics == 0 || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
2082
        if (dynamics)
2083
            delete dynamics;
peastman's avatar
peastman committed
2084
        dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), stepSize);
2085
        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
2086
2087
        prevStepSize = stepSize;
    }
2088
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
2089
    data.time += stepSize;
2090
    data.stepCount++;
2091
}
2092

2093
double ReferenceIntegrateVerletStepKernel::computeKineticEnergy(ContextImpl& context, const VerletIntegrator& integrator) {
2094
    return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
2095
2096
}

2097
2098
2099
2100
ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
}
2101

2102
void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
2103
    int numParticles = system.getNumParticles();
2104
    masses.resize(numParticles);
2105
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2106
        masses[i] = system.getParticleMass(i);
2107
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
2108
2109
}

2110
void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
2111
2112
2113
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
peastman's avatar
peastman committed
2114
2115
2116
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2117
2118
2119
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
2120
        if (dynamics)
2121
            delete dynamics;
2122
        dynamics = new ReferenceStochasticDynamics(
2123
                context.getSystem().getNumParticles(), 
peastman's avatar
peastman committed
2124
2125
2126
                stepSize, 
                friction, 
                temperature);
2127
        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
2128
2129
2130
2131
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
2132
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
2133
    data.time += stepSize;
2134
    data.stepCount++;
2135
2136
}

2137
double ReferenceIntegrateLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const LangevinIntegrator& integrator) {
2138
    return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
2139
2140
}

2141
2142
2143
2144
2145
ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
    if (dynamics)
        delete dynamics;
}

2146
void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
2147
    int numParticles = system.getNumParticles();
2148
    masses.resize(numParticles);
2149
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2150
        masses[i] = system.getParticleMass(i);
2151
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
2152
2153
}

2154
void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
2155
2156
2157
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
peastman's avatar
peastman committed
2158
2159
2160
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2161
2162
2163
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
2164
        if (dynamics)
2165
            delete dynamics;
2166
        dynamics = new ReferenceBrownianDynamics(
2167
                context.getSystem().getNumParticles(), 
peastman's avatar
peastman committed
2168
2169
2170
                stepSize, 
                friction, 
                temperature);
2171
        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
2172
2173
2174
2175
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
2176
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance());
2177
    data.time += stepSize;
2178
    data.stepCount++;
2179
2180
}

2181
double ReferenceIntegrateBrownianStepKernel::computeKineticEnergy(ContextImpl& context, const BrownianIntegrator& integrator) {
2182
    return computeShiftedKineticEnergy(context, masses, 0);
2183
2184
}

2185
2186
2187
2188
2189
2190
2191
ReferenceIntegrateVariableLangevinStepKernel::~ReferenceIntegrateVariableLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
}

void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& system, const VariableLangevinIntegrator& integrator) {
    int numParticles = system.getNumParticles();
2192
    masses.resize(numParticles);
2193
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2194
        masses[i] = system.getParticleMass(i);
2195
2196
2197
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
}

2198
double ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
2199
2200
2201
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double errorTol = integrator.getErrorTolerance();
peastman's avatar
peastman committed
2202
2203
2204
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2205
2206
2207
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || errorTol != prevErrorTol) {
        // Recreate the computation objects with the new parameters.

2208
        if (dynamics)
2209
            delete dynamics;
peastman's avatar
peastman committed
2210
        dynamics = new ReferenceVariableStochasticDynamics(context.getSystem().getNumParticles(), friction, temperature, errorTol);
2211
        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
2212
2213
2214
2215
        prevTemp = temperature;
        prevFriction = friction;
        prevErrorTol = errorTol;
    }
peastman's avatar
peastman committed
2216
    double maxStepSize = maxTime-data.time;
2217
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize, integrator.getConstraintTolerance());
2218
2219
2220
2221
    data.time += dynamics->getDeltaT();
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
2222
    return dynamics->getDeltaT();
2223
2224
}

2225
double ReferenceIntegrateVariableLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const VariableLangevinIntegrator& integrator) {
2226
    return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
2227
2228
}

2229
2230
2231
2232
2233
2234
2235
ReferenceIntegrateVariableVerletStepKernel::~ReferenceIntegrateVariableVerletStepKernel() {
    if (dynamics)
        delete dynamics;
}

void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system, const VariableVerletIntegrator& integrator) {
    int numParticles = system.getNumParticles();
2236
    masses.resize(numParticles);
2237
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2238
        masses[i] = system.getParticleMass(i);
2239
2240
}

2241
double ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
2242
    double errorTol = integrator.getErrorTolerance();
peastman's avatar
peastman committed
2243
2244
2245
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2246
    if (dynamics == 0 || errorTol != prevErrorTol) {
2247
2248
        // Recreate the computation objects with the new parameters.

2249
        if (dynamics)
2250
            delete dynamics;
peastman's avatar
peastman committed
2251
        dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), errorTol);
2252
        dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
2253
        prevErrorTol = errorTol;
2254
    }
peastman's avatar
peastman committed
2255
    double maxStepSize = maxTime-data.time;
2256
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize, integrator.getConstraintTolerance());
2257
    data.time += dynamics->getDeltaT();
2258
2259
2260
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
2261
    return dynamics->getDeltaT();
2262
2263
}

2264
double ReferenceIntegrateVariableVerletStepKernel::computeKineticEnergy(ContextImpl& context, const VariableVerletIntegrator& integrator) {
2265
    return computeShiftedKineticEnergy(context, masses, 0.5*integrator.getStepSize());
2266
2267
}

2268
2269
2270
2271
2272
2273
2274
2275
2276
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
2277
        masses[i] = system.getParticleMass(i);
2278
    perDofValues.resize(integrator.getNumPerDofVariables());
peastman's avatar
peastman committed
2279
2280
    for (auto& values : perDofValues)
        values.resize(numParticles);
2281
2282
2283
2284

    // Create the computation objects.

    dynamics = new ReferenceCustomDynamics(system.getNumParticles(), integrator);
2285
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
2286
2287
2288
}

void ReferenceIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid) {
peastman's avatar
peastman committed
2289
2290
2291
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2292
2293
2294
2295
2296
2297
2298
2299
2300
2301
    
    // 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.
    
2302
2303
    dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
    dynamics->update(context, context.getSystem().getNumParticles(), posData, velData, forceData, masses, globals, perDofValues, forcesAreValid, integrator.getConstraintTolerance());
2304
2305
2306
2307
2308
2309
2310
2311
2312
2313
    
    // 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++;
}

2314
double ReferenceIntegrateCustomStepKernel::computeKineticEnergy(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid) {
peastman's avatar
peastman committed
2315
2316
2317
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& velData = extractVelocities(context);
    vector<Vec3>& forceData = extractForces(context);
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328
2329
2330
    
    // 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);
}

2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350
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];
}

2351
2352
2353
2354
2355
ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
    if (thermostat)
        delete thermostat;
}

2356
void ReferenceApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
Peter Eastman's avatar
Peter Eastman committed
2357
    int numParticles = system.getNumParticles();
2358
    masses.resize(numParticles);
2359
    for (int i = 0; i < numParticles; ++i)
peastman's avatar
peastman committed
2360
        masses[i] = system.getParticleMass(i);
2361
    this->thermostat = new ReferenceAndersenThermostat();
2362
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
2363
    particleGroups = AndersenThermostatImpl::calcParticleGroups(system);
2364
2365
}

2366
void ReferenceApplyAndersenThermostatKernel::execute(ContextImpl& context) {
peastman's avatar
peastman committed
2367
    vector<Vec3>& velData = extractVelocities(context);
2368
    thermostat->applyThermostat(particleGroups, velData, masses,
peastman's avatar
peastman committed
2369
2370
2371
        context.getParameter(AndersenThermostat::Temperature()),
        context.getParameter(AndersenThermostat::CollisionFrequency()),
        context.getIntegrator().getStepSize());
2372
2373
}

2374
2375
2376
2377
2378
ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel() {
    if (barostat)
        delete barostat;
}

2379
void ReferenceApplyMonteCarloBarostatKernel::initialize(const System& system, const Force& barostat) {
2380
2381
}

2382
void ReferenceApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ) {
2383
2384
    if (barostat == NULL)
        barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), context.getMolecules());
peastman's avatar
peastman committed
2385
2386
    vector<Vec3>& posData = extractPositions(context);
    Vec3* boxVectors = extractBoxVectors(context);
2387
    barostat->applyBarostat(posData, boxVectors, scaleX, scaleY, scaleZ);
2388
2389
2390
}

void ReferenceApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
peastman's avatar
peastman committed
2391
    vector<Vec3>& posData = extractPositions(context);
2392
2393
2394
    barostat->restorePositions(posData);
}

2395
2396
void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
    frequency = force.getFrequency();
Peter Eastman's avatar
Peter Eastman committed
2397
    masses.resize(system.getNumParticles());
2398
    for (size_t i = 0; i < masses.size(); ++i)
Peter Eastman's avatar
Peter Eastman committed
2399
        masses[i] = system.getParticleMass(i);
2400
2401
}

2402
void ReferenceRemoveCMMotionKernel::execute(ContextImpl& context) {
2403
    if (data.stepCount%frequency != 0)
2404
        return;
peastman's avatar
peastman committed
2405
    vector<Vec3>& velData = extractVelocities(context);
2406
2407
2408
    
    // Calculate the center of mass momentum.
    
peastman's avatar
peastman committed
2409
2410
    double momentum[] = {0.0, 0.0, 0.0};
    double mass = 0.0;
2411
    for (size_t i = 0; i < masses.size(); ++i) {
peastman's avatar
peastman committed
2412
2413
2414
2415
        momentum[0] += masses[i]*velData[i][0];
        momentum[1] += masses[i]*velData[i][1];
        momentum[2] += masses[i]*velData[i][2];
        mass += masses[i];
2416
2417
    }
    
Peter Eastman's avatar
Peter Eastman committed
2418
    // Adjust the particle velocities.
2419
    
2420
2421
2422
    momentum[0] /= mass;
    momentum[1] /= mass;
    momentum[2] /= mass;
2423
    for (size_t i = 0; i < masses.size(); ++i) {
2424
2425
2426
2427
2428
        if (masses[i] != 0.0) {
            velData[i][0] -= momentum[0];
            velData[i][1] -= momentum[1];
            velData[i][2] -= momentum[2];
        }
2429
2430
    }
}