ReferenceKernels.cpp 82.7 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-2009 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 "gbsa/CpuObc.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
34
#include "gbsa/CpuGBVI.h"
35
#include "SimTKReference/ReferenceAndersenThermostat.h"
36
37
#include "SimTKReference/ReferenceAngleBondIxn.h"
#include "SimTKReference/ReferenceBondForce.h"
38
#include "SimTKReference/ReferenceBrownianDynamics.h"
39
#include "SimTKReference/ReferenceCCMAAlgorithm.h"
40
#include "SimTKReference/ReferenceCMAPTorsionIxn.h"
41
#include "SimTKReference/ReferenceCustomAngleIxn.h"
42
#include "SimTKReference/ReferenceCustomBondIxn.h"
43
#include "SimTKReference/ReferenceCustomExternalIxn.h"
44
#include "SimTKReference/ReferenceCustomGBIxn.h"
45
#include "SimTKReference/ReferenceCustomHbondIxn.h"
46
#include "SimTKReference/ReferenceCustomNonbondedIxn.h"
47
#include "SimTKReference/ReferenceCustomTorsionIxn.h"
48
49
50
#include "SimTKReference/ReferenceHarmonicBondIxn.h"
#include "SimTKReference/ReferenceLJCoulomb14.h"
#include "SimTKReference/ReferenceLJCoulombIxn.h"
51
#include "SimTKReference/ReferenceMonteCarloBarostat.h"
52
53
#include "SimTKReference/ReferenceProperDihedralBond.h"
#include "SimTKReference/ReferenceRbDihedralBond.h"
54
#include "SimTKReference/ReferenceStochasticDynamics.h"
55
56
#include "SimTKReference/ReferenceVariableStochasticDynamics.h"
#include "SimTKReference/ReferenceVariableVerletDynamics.h"
57
#include "SimTKReference/ReferenceVerletDynamics.h"
58
#include "openmm/CMMotionRemover.h"
59
#include "openmm/Context.h"
60
#include "openmm/System.h"
61
#include "openmm/internal/ContextImpl.h"
62
#include "openmm/internal/CustomHbondForceImpl.h"
63
#include "openmm/internal/CMAPTorsionForceImpl.h"
64
#include "openmm/internal/NonbondedForceImpl.h"
65
#include "openmm/Integrator.h"
66
#include "openmm/OpenMMException.h"
67
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
68
#include "lepton/CustomFunction.h"
69
70
#include "lepton/Parser.h"
#include "lepton/ParsedExpression.h"
71
#include <cmath>
72
#include <limits>
73
74
75
76

using namespace OpenMM;
using namespace std;

77
static int** allocateIntArray(int length, int width) {
78
79
80
81
82
83
    int** array = new int*[length];
    for (int i = 0; i < length; ++i)
        array[i] = new int[width];
    return array;
}

84
static RealOpenMM** allocateRealArray(int length, int width) {
85
86
87
88
89
90
    RealOpenMM** array = new RealOpenMM*[length];
    for (int i = 0; i < length; ++i)
        array[i] = new RealOpenMM[width];
    return array;
}

91
static int** copyToArray(const vector<vector<int> > vec) {
92
93
94
    if (vec.size() == 0)
        return new int*[0];
    int** array = allocateIntArray(vec.size(), vec[0].size());
95
96
    for (size_t i = 0; i < vec.size(); ++i)
        for (size_t j = 0; j < vec[i].size(); ++j)
97
98
99
100
            array[i][j] = vec[i][j];
    return array;
}

101
static RealOpenMM** copyToArray(const vector<vector<double> > vec) {
102
103
104
    if (vec.size() == 0)
        return new RealOpenMM*[0];
    RealOpenMM** array = allocateRealArray(vec.size(), vec[0].size());
105
106
107
    for (size_t i = 0; i < vec.size(); ++i)
        for (size_t j = 0; j < vec[i].size(); ++j)
            array[i][j] = static_cast<RealOpenMM>(vec[i][j]);
108
109
110
    return array;
}

111
static void disposeIntArray(int** array, int size) {
112
113
114
115
116
117
118
    if (array) {
        for (int i = 0; i < size; ++i)
            delete[] array[i];
        delete[] array;
    }
}

119
static void disposeRealArray(RealOpenMM** array, int size) {
120
121
122
123
124
125
126
    if (array) {
        for (int i = 0; i < size; ++i)
            delete[] array[i];
        delete[] array;
    }
}

127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
static RealOpenMM** extractPositions(ContextImpl& context) {
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
    return (RealOpenMM**) data->positions;
}

static RealOpenMM** extractVelocities(ContextImpl& context) {
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
    return (RealOpenMM**) data->velocities;
}

static RealOpenMM** extractForces(ContextImpl& context) {
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
    return (RealOpenMM**) data->forces;
}

142
143
144
145
146
static RealOpenMM* extractBoxSize(ContextImpl& context) {
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
    return (RealOpenMM*) data->periodicBoxSize;
}

147
static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorithm::AngleInfo>& angles) {
148
149
150
151
152
153
154
    for (int i = 0; i < system.getNumForces(); i++) {
        const HarmonicAngleForce* force = dynamic_cast<const HarmonicAngleForce*>(&system.getForce(i));
        if (force != NULL) {
            for (int j = 0; j < force->getNumAngles(); j++) {
                int atom1, atom2, atom3;
                double angle, k;
                force->getAngleParameters(j, atom1, atom2, atom3, angle, k);
155
                angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, (RealOpenMM)angle));
156
157
158
159
160
            }
        }
    }
}

161
void ReferenceCalcForcesAndEnergyKernel::initialize(const System& system) {
162
163
}

164
void ReferenceCalcForcesAndEnergyKernel::beginForceComputation(ContextImpl& context) {
165
166
167
168
169
170
171
    int numParticles = context.getSystem().getNumParticles();
    RealOpenMM** forceData = extractForces(context);
    for (int i = 0; i < numParticles; ++i) {
        forceData[i][0] = (RealOpenMM) 0.0;
        forceData[i][1] = (RealOpenMM) 0.0;
        forceData[i][2] = (RealOpenMM) 0.0;
    }
172
173
}

174
175
176
177
178
179
180
181
182
183
void ReferenceCalcForcesAndEnergyKernel::finishForceComputation(ContextImpl& context) {
}

void ReferenceCalcForcesAndEnergyKernel::beginEnergyComputation(ContextImpl& context) {
}

double ReferenceCalcForcesAndEnergyKernel::finishEnergyComputation(ContextImpl& context) {
    return 0.0;
}

184
void ReferenceUpdateStateDataKernel::initialize(const System& system) {
185
186
}

187
double ReferenceUpdateStateDataKernel::getTime(const ContextImpl& context) const {
188
189
190
    return data.time;
}

191
void ReferenceUpdateStateDataKernel::setTime(ContextImpl& context, double time) {
192
193
194
    data.time = time;
}

195
196
197
198
199
200
201
202
203
204
205
206
void ReferenceUpdateStateDataKernel::getPositions(ContextImpl& context, std::vector<Vec3>& positions) {
    int numParticles = context.getSystem().getNumParticles();
    RealOpenMM** posData = extractPositions(context);
    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();
    RealOpenMM** posData = extractPositions(context);
    for (int i = 0; i < numParticles; ++i) {
207
208
209
        posData[i][0] = (RealOpenMM) positions[i][0];
        posData[i][1] = (RealOpenMM) positions[i][1];
        posData[i][2] = (RealOpenMM) positions[i][2];
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
    }
}

void ReferenceUpdateStateDataKernel::getVelocities(ContextImpl& context, std::vector<Vec3>& velocities) {
    int numParticles = context.getSystem().getNumParticles();
    RealOpenMM** velData = extractVelocities(context);
    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();
    RealOpenMM** velData = extractVelocities(context);
    for (int i = 0; i < numParticles; ++i) {
225
226
227
        velData[i][0] = (RealOpenMM) velocities[i][0];
        velData[i][1] = (RealOpenMM) velocities[i][1];
        velData[i][2] = (RealOpenMM) velocities[i][2];
228
229
230
231
232
233
234
235
236
237
238
    }
}

void ReferenceUpdateStateDataKernel::getForces(ContextImpl& context, std::vector<Vec3>& forces) {
    int numParticles = context.getSystem().getNumParticles();
    RealOpenMM** forceData = extractForces(context);
    forces.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        forces[i] = Vec3(forceData[i][0], forceData[i][1], forceData[i][2]);
}

239
240
241
242
243
244
245
246
247
248
249
250
251
252
void ReferenceUpdateStateDataKernel::getPeriodicBoxVectors(ContextImpl& context, Vec3& a, Vec3& b, Vec3& c) const {
    RealOpenMM* box = extractBoxSize(context);
    a = Vec3(box[0], 0, 0);
    b = Vec3(0, box[1], 0);
    c = Vec3(0, 0, box[2]);
}

void ReferenceUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c) const {
    RealOpenMM* box = extractBoxSize(context);
    box[0] = (RealOpenMM) a[0];
    box[1] = (RealOpenMM) b[1];
    box[2] = (RealOpenMM) c[2];
}

253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
void ReferenceApplyConstraintsKernel::initialize(const System& system) {
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
    inverseMasses = new RealOpenMM[numParticles];
    for (int i = 0; i < numParticles; ++i) {
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
        inverseMasses[i] = 1.0/masses[i];
    }
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
    constraintDistances = new RealOpenMM[numConstraints];
    for (int i = 0; i < numConstraints; ++i) {
        int particle1, particle2;
        double distance;
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
    }
}

ReferenceApplyConstraintsKernel::~ReferenceApplyConstraintsKernel() {
    if (constraints)
        delete constraints;
    if (masses)
        delete[] masses;
    if (inverseMasses)
        delete[] inverseMasses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
    if (constraintDistances)
        delete[] constraintDistances;
}

void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
    if (constraints == NULL) {
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, tol);
    }
    RealOpenMM** positions = extractPositions(context);
    constraints->setTolerance(tol);
    constraints->apply(data.numParticles, positions, positions, inverseMasses);
}

298
ReferenceCalcHarmonicBondForceKernel::~ReferenceCalcHarmonicBondForceKernel() {
299
300
301
302
    disposeIntArray(bondIndexArray, numBonds);
    disposeRealArray(bondParamArray, numBonds);
}

303
void ReferenceCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
304
305
306
307
    numBonds = force.getNumBonds();
    bondIndexArray = allocateIntArray(numBonds, 2);
    bondParamArray = allocateRealArray(numBonds, 2);
    for (int i = 0; i < force.getNumBonds(); ++i) {
Peter Eastman's avatar
Peter Eastman committed
308
        int particle1, particle2;
309
        double length, k;
Peter Eastman's avatar
Peter Eastman committed
310
311
312
        force.getBondParameters(i, particle1, particle2, length, k);
        bondIndexArray[i][0] = particle1;
        bondIndexArray[i][1] = particle2;
313
314
        bondParamArray[i][0] = (RealOpenMM) length;
        bondParamArray[i][1] = (RealOpenMM) k;
315
    }
316
317
}

318
void ReferenceCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
319
320
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
321
322
    ReferenceBondForce refBondForce;
    ReferenceHarmonicBondIxn harmonicBond;
323
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, harmonicBond);
324
325
}

326
double ReferenceCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
327
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
328
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
329
330
331
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceHarmonicBondIxn harmonicBond;
332
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, &energy, harmonicBond);
Peter Eastman's avatar
Peter Eastman committed
333
    disposeRealArray(forceData, context.getSystem().getNumParticles());
334
335
336
    return energy;
}

337
338
339
340
341
342
343
344
345
346
347
ReferenceCalcCustomBondForceKernel::~ReferenceCalcCustomBondForceKernel() {
    disposeIntArray(bondIndexArray, numBonds);
    disposeRealArray(bondParamArray, numBonds);
}

void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const CustomBondForce& force) {
    numBonds = force.getNumBonds();
    int numParameters = force.getNumPerBondParameters();

    // Build the arrays.

348
    bondIndexArray = allocateIntArray(numBonds, 2);
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
    bondParamArray = allocateRealArray(numBonds, numParameters);
    vector<double> params;
    for (int i = 0; i < force.getNumBonds(); ++i) {
        int particle1, particle2;
        force.getBondParameters(i, particle1, particle2, params);
        bondIndexArray[i][0] = particle1;
        bondIndexArray[i][1] = particle2;
        for (int j = 0; j < numParameters; j++)
            bondParamArray[i][j] = (RealOpenMM) params[j];
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
    energyExpression = expression.createProgram();
    forceExpression = expression.differentiate("r").optimize().createProgram();
    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));
}

void ReferenceCalcCustomBondForceKernel::executeForces(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceBondForce refBondForce;
    ReferenceCustomBondIxn harmonicBond(energyExpression, forceExpression, parameterNames, globalParameters);
379
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, harmonicBond);
380
381
382
383
384
385
386
387
388
389
390
}

double ReferenceCalcCustomBondForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
    RealOpenMM energy = 0;
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceBondForce refBondForce;
    ReferenceCustomBondIxn harmonicBond(energyExpression, forceExpression, parameterNames, globalParameters);
391
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, &energy, harmonicBond);
392
393
394
395
    disposeRealArray(forceData, context.getSystem().getNumParticles());
    return energy;
}

396
397
398
399
400
401
402
403
404
ReferenceCalcHarmonicAngleForceKernel::~ReferenceCalcHarmonicAngleForceKernel() {
    disposeIntArray(angleIndexArray, numAngles);
    disposeRealArray(angleParamArray, numAngles);
}

void ReferenceCalcHarmonicAngleForceKernel::initialize(const System& system, const HarmonicAngleForce& force) {
    numAngles = force.getNumAngles();
    angleIndexArray = allocateIntArray(numAngles, 3);
    angleParamArray = allocateRealArray(numAngles, 2);
405
    for (int i = 0; i < force.getNumAngles(); ++i) {
Peter Eastman's avatar
Peter Eastman committed
406
        int particle1, particle2, particle3;
407
        double angle, k;
Peter Eastman's avatar
Peter Eastman committed
408
409
410
411
        force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
        angleIndexArray[i][0] = particle1;
        angleIndexArray[i][1] = particle2;
        angleIndexArray[i][2] = particle3;
412
413
        angleParamArray[i][0] = (RealOpenMM) angle;
        angleParamArray[i][1] = (RealOpenMM) k;
414
    }
415
416
}

417
void ReferenceCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
418
419
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
420
421
    ReferenceBondForce refBondForce;
    ReferenceAngleBondIxn angleBond;
422
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, angleBond);
423
424
}

425
double ReferenceCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
426
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
427
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
428
429
430
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceAngleBondIxn angleBond;
431
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, &energy, angleBond);
Peter Eastman's avatar
Peter Eastman committed
432
    disposeRealArray(forceData, context.getSystem().getNumParticles());
433
434
435
    return energy;
}

436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
ReferenceCalcCustomAngleForceKernel::~ReferenceCalcCustomAngleForceKernel() {
    disposeIntArray(angleIndexArray, numAngles);
    disposeRealArray(angleParamArray, numAngles);
}

void ReferenceCalcCustomAngleForceKernel::initialize(const System& system, const CustomAngleForce& force) {
    numAngles = force.getNumAngles();
    int numParameters = force.getNumPerAngleParameters();

    // Build the arrays.

    angleIndexArray = allocateIntArray(numAngles, 3);
    angleParamArray = allocateRealArray(numAngles, numParameters);
    vector<double> params;
    for (int i = 0; i < force.getNumAngles(); ++i) {
        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++)
            angleParamArray[i][j] = (RealOpenMM) params[j];
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
    energyExpression = expression.createProgram();
    forceExpression = expression.differentiate("theta").optimize().createProgram();
    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));
}

void ReferenceCalcCustomAngleForceKernel::executeForces(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceBondForce refBondForce;
478
    ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters);
479
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, customAngle);
480
481
482
483
484
485
486
487
488
489
}

double ReferenceCalcCustomAngleForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
    RealOpenMM energy = 0;
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceBondForce refBondForce;
490
    ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters);
491
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, &energy, customAngle);
492
493
494
495
    disposeRealArray(forceData, context.getSystem().getNumParticles());
    return energy;
}

496
497
498
499
500
501
502
503
504
505
ReferenceCalcPeriodicTorsionForceKernel::~ReferenceCalcPeriodicTorsionForceKernel() {
    disposeIntArray(torsionIndexArray, numTorsions);
    disposeRealArray(torsionParamArray, numTorsions);
}

void ReferenceCalcPeriodicTorsionForceKernel::initialize(const System& system, const PeriodicTorsionForce& force) {
    numTorsions = force.getNumTorsions();
    torsionIndexArray = allocateIntArray(numTorsions, 4);
    torsionParamArray = allocateRealArray(numTorsions, 3);
    for (int i = 0; i < force.getNumTorsions(); ++i) {
Peter Eastman's avatar
Peter Eastman committed
506
        int particle1, particle2, particle3, particle4, periodicity;
507
        double phase, k;
Peter Eastman's avatar
Peter Eastman committed
508
509
510
511
512
        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;
513
514
515
        torsionParamArray[i][0] = (RealOpenMM) k;
        torsionParamArray[i][1] = (RealOpenMM) phase;
        torsionParamArray[i][2] = (RealOpenMM) periodicity;
516
    }
517
518
}

519
void ReferenceCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
520
521
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
522
523
    ReferenceBondForce refBondForce;
    ReferenceProperDihedralBond periodicTorsionBond;
524
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, periodicTorsionBond);
525
526
}

527
double ReferenceCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
528
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
529
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
530
531
532
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceProperDihedralBond periodicTorsionBond;
533
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, &energy, periodicTorsionBond);
Peter Eastman's avatar
Peter Eastman committed
534
    disposeRealArray(forceData, context.getSystem().getNumParticles());
535
536
537
538
539
540
541
542
543
544
545
546
547
    return energy;
}

ReferenceCalcRBTorsionForceKernel::~ReferenceCalcRBTorsionForceKernel() {
    disposeIntArray(torsionIndexArray, numTorsions);
    disposeRealArray(torsionParamArray, numTorsions);
}

void ReferenceCalcRBTorsionForceKernel::initialize(const System& system, const RBTorsionForce& force) {
    numTorsions = force.getNumTorsions();
    torsionIndexArray = allocateIntArray(numTorsions, 4);
    torsionParamArray = allocateRealArray(numTorsions, 6);
    for (int i = 0; i < force.getNumTorsions(); ++i) {
Peter Eastman's avatar
Peter Eastman committed
548
        int particle1, particle2, particle3, particle4;
549
        double c0, c1, c2, c3, c4, c5;
Peter Eastman's avatar
Peter Eastman committed
550
551
552
553
554
        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;
555
556
557
558
559
560
        torsionParamArray[i][0] = (RealOpenMM) c0;
        torsionParamArray[i][1] = (RealOpenMM) c1;
        torsionParamArray[i][2] = (RealOpenMM) c2;
        torsionParamArray[i][3] = (RealOpenMM) c3;
        torsionParamArray[i][4] = (RealOpenMM) c4;
        torsionParamArray[i][5] = (RealOpenMM) c5;
561
    }
562
563
}

564
void ReferenceCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
565
566
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
567
568
    ReferenceBondForce refBondForce;
    ReferenceRbDihedralBond rbTorsionBond;
569
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, rbTorsionBond);
570
571
}

572
double ReferenceCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
573
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
574
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
575
576
577
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceRbDihedralBond rbTorsionBond;
578
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, &energy, rbTorsionBond);
Peter Eastman's avatar
Peter Eastman committed
579
    disposeRealArray(forceData, context.getSystem().getNumParticles());
580
581
582
    return energy;
}

583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
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]);
    }
}

void ReferenceCalcCMAPTorsionForceKernel::executeForces(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
    ReferenceCMAPTorsionIxn torsion(coeff, torsionMaps, torsionIndices);
    torsion.calculateIxn(posData, forceData, 0);
}

double ReferenceCalcCMAPTorsionForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
    RealOpenMM totalEnergy = 0;
    ReferenceCMAPTorsionIxn torsion(coeff, torsionMaps, torsionIndices);
    torsion.calculateIxn(posData, forceData, &totalEnergy);
    return totalEnergy;
}

625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
ReferenceCalcCustomTorsionForceKernel::~ReferenceCalcCustomTorsionForceKernel() {
    disposeIntArray(torsionIndexArray, numTorsions);
    disposeRealArray(torsionParamArray, numTorsions);
}

void ReferenceCalcCustomTorsionForceKernel::initialize(const System& system, const CustomTorsionForce& force) {
    numTorsions = force.getNumTorsions();
    int numParameters = force.getNumPerTorsionParameters();

    // Build the arrays.

    torsionIndexArray = allocateIntArray(numTorsions, 4);
    torsionParamArray = allocateRealArray(numTorsions, numParameters);
    vector<double> params;
    for (int i = 0; i < force.getNumTorsions(); ++i) {
        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++)
            torsionParamArray[i][j] = (RealOpenMM) params[j];
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
    energyExpression = expression.createProgram();
    forceExpression = expression.differentiate("theta").optimize().createProgram();
    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));
}

void ReferenceCalcCustomTorsionForceKernel::executeForces(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceBondForce refBondForce;
    ReferenceCustomTorsionIxn customTorsion(energyExpression, forceExpression, parameterNames, globalParameters);
669
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, customTorsion);
670
671
672
673
674
675
676
677
678
679
680
}

double ReferenceCalcCustomTorsionForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
    RealOpenMM energy = 0;
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceBondForce refBondForce;
    ReferenceCustomTorsionIxn customTorsion(energyExpression, forceExpression, parameterNames, globalParameters);
681
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, &energy, customTorsion);
682
683
684
685
    disposeRealArray(forceData, context.getSystem().getNumParticles());
    return energy;
}

686
ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
Peter Eastman's avatar
Peter Eastman committed
687
688
    disposeRealArray(particleParamArray, numParticles);
    disposeIntArray(exclusionArray, numParticles);
689
690
691
692
693
694
    disposeIntArray(bonded14IndexArray, num14);
    disposeRealArray(bonded14ParamArray, num14);
    if (neighborList != NULL)
        delete neighborList;
}

695
696
697
698
void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {

    // Identify which exceptions are 1-4 interactions.

Peter Eastman's avatar
Peter Eastman committed
699
    numParticles = force.getNumParticles();
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
    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();
715
716
    bonded14IndexArray = allocateIntArray(num14, 2);
    bonded14ParamArray = allocateRealArray(num14, 3);
Peter Eastman's avatar
Peter Eastman committed
717
718
    particleParamArray = allocateRealArray(numParticles, 3);
    for (int i = 0; i < numParticles; ++i) {
719
        double charge, radius, depth;
Peter Eastman's avatar
Peter Eastman committed
720
721
722
        force.getParticleParameters(i, charge, radius, depth);
        particleParamArray[i][0] = static_cast<RealOpenMM>(0.5*radius);
        particleParamArray[i][1] = static_cast<RealOpenMM>(2.0*sqrt(depth));
723
        particleParamArray[i][2] = static_cast<RealOpenMM>(charge);
724
    }
725
    this->exclusions = exclusions;
Peter Eastman's avatar
Peter Eastman committed
726
727
    exclusionArray = new int*[numParticles];
    for (int i = 0; i < numParticles; ++i) {
728
729
730
731
732
733
734
        exclusionArray[i] = new int[exclusions[i].size()+1];
        exclusionArray[i][0] = exclusions[i].size();
        int index = 0;
        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
            exclusionArray[i][++index] = *iter;
    }
    for (int i = 0; i < num14; ++i) {
Peter Eastman's avatar
Peter Eastman committed
735
        int particle1, particle2;
736
        double charge, radius, depth;
737
        force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth);
Peter Eastman's avatar
Peter Eastman committed
738
739
        bonded14IndexArray[i][0] = particle1;
        bonded14IndexArray[i][1] = particle2;
740
741
        bonded14ParamArray[i][0] = static_cast<RealOpenMM>(radius);
        bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
742
        bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge);
743
    }
744
    nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
745
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
746
747
748
749
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();
750
751
752
    if (nonbondedMethod == Ewald) {
        double alpha;
        NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmax[0], kmax[1], kmax[2]);
753
        ewaldAlpha = (RealOpenMM) alpha;
754
755
756
757
    }
    else if (nonbondedMethod == PME) {
        double alpha;
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2]);
758
        ewaldAlpha = (RealOpenMM) alpha;
759
    }
760
    rfDielectric = (RealOpenMM)force.getReactionFieldDielectric();
761
762
763
764
    if (force.getUseDispersionCorrection())
        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
    else
        dispersionCoefficient = 0.0;
765
766
}

767
void ReferenceCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
768
769
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
770
    ReferenceLJCoulombIxn clj;
771
    bool periodic = (nonbondedMethod == CutoffPeriodic);
772
    bool ewald  = (nonbondedMethod == Ewald);
773
    bool pme  = (nonbondedMethod == PME);
774
    if (nonbondedMethod != NoCutoff) {
775
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
776
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
777
    }
778
    if (periodic || ewald || pme)
779
        clj.setPeriodic(extractBoxSize(context));
780
781
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
782
    if (pme)
783
        clj.setUsePME(ewaldAlpha, gridSize);
784
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
785
    ReferenceBondForce refBondForce;
786
    ReferenceLJCoulomb14 nonbonded14;
787
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, nonbonded14);
788
789
}

790
double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
791
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
792
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
793
794
    RealOpenMM energy = 0;
    ReferenceLJCoulombIxn clj;
795
    bool periodic = (nonbondedMethod == CutoffPeriodic);
796
    bool ewald  = (nonbondedMethod == Ewald);
797
    bool pme  = (nonbondedMethod == PME);
798
    if (nonbondedMethod != NoCutoff) {
799
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
800
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
801
    }
802
    if (periodic || ewald || pme)
803
        clj.setPeriodic(extractBoxSize(context));
804
805
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
806
    if (pme)
807
        clj.setUsePME(ewaldAlpha, gridSize);
808
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
809
    ReferenceBondForce refBondForce;
810
    ReferenceLJCoulomb14 nonbonded14;
811
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, &energy, nonbonded14);
Peter Eastman's avatar
Peter Eastman committed
812
    disposeRealArray(forceData, numParticles);
813
814
815
816
    if (periodic || ewald || pme) {
        RealOpenMM* boxSize = extractBoxSize(context);
        energy += dispersionCoefficient/(boxSize[0]*boxSize[1]*boxSize[2]);
    }
817
818
819
    return energy;
}

820
class ReferenceTabulatedFunction : public Lepton::CustomFunction {
821
public:
822
    ReferenceTabulatedFunction(double min, double max, const vector<double>& values, bool interpolating) :
823
824
825
826
827
828
829
830
831
832
            min(min), max(max), values(values), interpolating(interpolating) {
    }
    int getNumArguments() const {
        return 1;
    }
    /**
     * Given the function argument, find the local spline coefficients.
     */
    void findCoefficients(double& x, double* coeff) const {
        int length = values.size();
833
        double scale = (length-1)/(max-min);
834
        int index = (int) std::floor((x-min)*scale);
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
        double points[4];
        points[0] = (index == 0 ? 2*values[0]-values[1] : values[index-1]);
        points[1] = values[index];
        points[2] = (index > length-2 ? values[length-1] : values[index+1]);
        points[3] = (index > length-3 ? 2*values[length-1]-values[length-2] : values[index+2]);
        if (interpolating) {
            coeff[0] = points[1];
            coeff[1] = 0.5*(-points[0]+points[2]);
            coeff[2] = 0.5*(2.0*points[0]-5.0*points[1]+4.0*points[2]-points[3]);
            coeff[3] = 0.5*(-points[0]+3.0*points[1]-3.0*points[2]+points[3]);
        }
        else {
            coeff[0] = (points[0]+4.0*points[1]+points[2])/6.0;
            coeff[1] = (-3.0*points[0]+3.0*points[2])/6.0;
            coeff[2] = (3.0*points[0]-6.0*points[1]+3.0*points[2])/6.0;
            coeff[3] = (-points[0]+3.0*points[1]-3.0*points[2]+points[3])/6.0;
        }
852
        x = (x-min)*scale-index;
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
    }
    double evaluate(const double* arguments) const {
        double x = arguments[0];
        if (x < min || x > max)
            return 0.0;
        double coeff[4];
        findCoefficients(x, coeff);
        return coeff[0]+x*(coeff[1]+x*(coeff[2]+x*coeff[3]));
    }
    double evaluateDerivative(const double* arguments, const int* derivOrder) const {
        double x = arguments[0];
        if (x < min || x > max)
            return 0.0;
        double coeff[4];
        findCoefficients(x, coeff);
        double scale = (values.size()-1)/(max-min);
869
        return scale*(coeff[1]+x*(2.0*coeff[2]+x*3.0*coeff[3])); // We assume a first derivative, because that's the only order ever used by CustomNonbondedForce.
870
871
    }
    CustomFunction* clone() const {
872
        return new ReferenceTabulatedFunction(min, max, values, interpolating);
873
874
875
876
877
878
    }
    double min, max;
    vector<double> values;
    bool interpolating;
};

879
880
881
882
883
884
885
886
887
ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
    disposeRealArray(particleParamArray, numParticles);
    disposeIntArray(exclusionArray, numParticles);
    if (neighborList != NULL)
        delete neighborList;
}

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

888
    // Record the exclusions.
889
890
891

    numParticles = force.getNumParticles();
    exclusions.resize(numParticles);
892
    for (int i = 0; i < force.getNumExclusions(); i++) {
893
        int particle1, particle2;
894
        force.getExclusionParticles(i, particle1, particle2);
895
896
897
898
899
900
        exclusions[particle1].insert(particle2);
        exclusions[particle2].insert(particle1);
    }

    // Build the arrays.

901
    int numParameters = force.getNumPerParticleParameters();
902
903
    particleParamArray = allocateRealArray(numParticles, numParameters);
    for (int i = 0; i < numParticles; ++i) {
904
        vector<double> parameters;
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
        force.getParticleParameters(i, parameters);
        for (int j = 0; j < numParameters; j++)
            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
    }
    exclusionArray = new int*[numParticles];
    for (int i = 0; i < numParticles; ++i) {
        exclusionArray[i] = new int[exclusions[i].size()+1];
        exclusionArray[i][0] = exclusions[i].size();
        int index = 0;
        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
            exclusionArray[i][++index] = *iter;
    }
    nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();

924
925
926
927
928
929
930
931
932
    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
    for (int i = 0; i < force.getNumFunctions(); i++) {
        string name;
        vector<double> values;
        double min, max;
        bool interpolating;
        force.getFunctionParameters(i, name, values, min, max, interpolating);
933
        functions[name] = new ReferenceTabulatedFunction(min, max, values, interpolating);
934
935
    }

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

938
    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
939
940
    energyExpression = expression.createProgram();
    forceExpression = expression.differentiate("r").optimize().createProgram();
941
942
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerParticleParameterName(i));
943
944
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
945
946
947
948
949

    // Delete the custom functions.

    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
        delete iter->second;
950
951
952
}

void ReferenceCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
953
954
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
955
    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames);
956
957
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
958
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
959
960
961
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    if (periodic)
962
        ixn.setPeriodic(extractBoxSize(context));
963
    map<string, double> globalParameters;
964
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
965
966
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, 0);
967
968
969
}

double ReferenceCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
970
    RealOpenMM** posData = extractPositions(context);
971
972
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
    RealOpenMM energy = 0;
973
    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames);
974
975
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
976
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
977
978
979
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    if (periodic)
980
        ixn.setPeriodic(extractBoxSize(context));
981
    map<string, double> globalParameters;
982
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
983
984
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, &energy);
985
986
987
988
    disposeRealArray(forceData, numParticles);
    return energy;
}

989
ReferenceCalcGBSAOBCForceKernel::~ReferenceCalcGBSAOBCForceKernel() {
990
    if (obc) {
991
        // delete obc->getObcParameters();
992
993
994
995
        delete obc;
    }
}

996
void ReferenceCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
Peter Eastman's avatar
Peter Eastman committed
997
998
999
1000
1001
    int numParticles = system.getNumParticles();
    charges.resize(numParticles);
    vector<RealOpenMM> atomicRadii(numParticles);
    vector<RealOpenMM> scaleFactors(numParticles);
    for (int i = 0; i < numParticles; ++i) {
1002
        double charge, radius, scalingFactor;
Peter Eastman's avatar
Peter Eastman committed
1003
        force.getParticleParameters(i, charge, radius, scalingFactor);
1004
1005
1006
        charges[i] = static_cast<RealOpenMM>(charge);
        atomicRadii[i] = static_cast<RealOpenMM>(radius);
        scaleFactors[i] = static_cast<RealOpenMM>(scalingFactor);
1007
    }
1008
    ObcParameters* obcParameters = new ObcParameters(numParticles, ObcParameters::ObcTypeII);
1009
    obcParameters->setAtomicRadii(atomicRadii);
1010
    obcParameters->setScaledRadiusFactors(scaleFactors);
1011
1012
    obcParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
    obcParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
1013
1014
    if (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff)
        obcParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
1015
    isPeriodic = (force.getNonbondedMethod() == GBSAOBCForce::CutoffPeriodic);
1016
1017
    obc = new CpuObc(obcParameters);
    obc->setIncludeAceApproximation(true);
1018
1019
}

1020
void ReferenceCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
1021
1022
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
1023
1024
    if (isPeriodic)
        obc->getObcParameters()->setPeriodic(extractBoxSize(context));
1025
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
1026
1027
}

1028
double ReferenceCalcGBSAOBCForceKernel::executeEnergy(ContextImpl& context) {
1029
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
1030
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
1031
1032
    if (isPeriodic)
        obc->getObcParameters()->setPeriodic(extractBoxSize(context));
1033
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
Peter Eastman's avatar
Peter Eastman committed
1034
    disposeRealArray(forceData, context.getSystem().getNumParticles());
1035
    return obc->getEnergy();
1036
1037
}

Mark Friedrichs's avatar
Mark Friedrichs committed
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
ReferenceCalcGBVIForceKernel::~ReferenceCalcGBVIForceKernel() {
    if (gbvi) {
        delete gbvi;
    }
}

void ReferenceCalcGBVIForceKernel::initialize(const System& system, const GBVIForce& force, const std::vector<double> & inputScaledRadii ) {
    int numParticles = system.getNumParticles();
    charges.resize(numParticles);
    vector<RealOpenMM> atomicRadii(numParticles);
    vector<RealOpenMM> scaledRadii(numParticles);
    vector<RealOpenMM> gammas(numParticles);
    for (int i = 0; i < numParticles; ++i) {
        double charge, radius, gamma;
        force.getParticleParameters(i, charge, radius, gamma);
        charges[i]       = static_cast<RealOpenMM>(charge);
        atomicRadii[i]   = static_cast<RealOpenMM>(radius);
        gammas[i]        = static_cast<RealOpenMM>(gamma);
        scaledRadii[i]   = static_cast<RealOpenMM>(inputScaledRadii[i]);
    }
    GBVIParameters * gBVIParameters = new GBVIParameters(numParticles);
    gBVIParameters->setAtomicRadii(atomicRadii);
    gBVIParameters->setGammaParameters(gammas);
    gBVIParameters->setScaledRadii(scaledRadii);
1062
1063
1064
1065
    gBVIParameters->setSolventDielectric(static_cast<RealOpenMM>(force.getSolventDielectric()));
    gBVIParameters->setSoluteDielectric(static_cast<RealOpenMM>(force.getSoluteDielectric()));
    if (force.getNonbondedMethod() != GBVIForce::NoCutoff)
        gBVIParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
1066
    isPeriodic = (force.getNonbondedMethod() == GBVIForce::CutoffPeriodic);
Mark Friedrichs's avatar
Mark Friedrichs committed
1067
1068
1069
    gbvi = new CpuGBVI(gBVIParameters);
}

1070
void ReferenceCalcGBVIForceKernel::executeForces(ContextImpl& context) {
Mark Friedrichs's avatar
Mark Friedrichs committed
1071

1072
1073
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
Mark Friedrichs's avatar
Mark Friedrichs committed
1074
    RealOpenMM* bornRadii  = new RealOpenMM[context.getSystem().getNumParticles()];
1075
1076
    if (isPeriodic)
        gbvi->getGBVIParameters()->setPeriodic(extractBoxSize(context));
Mark Friedrichs's avatar
Mark Friedrichs committed
1077
1078
1079
1080
1081
    gbvi->computeBornRadii(posData, bornRadii, NULL ); 
    gbvi->computeBornForces(bornRadii, posData, &charges[0], forceData);
    delete[] bornRadii;
}

1082
double ReferenceCalcGBVIForceKernel::executeEnergy(ContextImpl& context) {
1083
    RealOpenMM** posData = extractPositions(context);
Mark Friedrichs's avatar
Mark Friedrichs committed
1084
    RealOpenMM* bornRadii = new RealOpenMM[context.getSystem().getNumParticles()];
1085
1086
    if (isPeriodic)
        gbvi->getGBVIParameters()->setPeriodic(extractBoxSize(context));
1087
    gbvi->computeBornRadii(posData, bornRadii, NULL );
Mark Friedrichs's avatar
Mark Friedrichs committed
1088
1089
1090
1091
1092
    RealOpenMM energy     = gbvi->computeBornEnergy(bornRadii ,posData, &charges[0]);
    delete[] bornRadii;
    return static_cast<double>(energy);
}

1093
1094
1095
1096
1097
1098
1099
ReferenceCalcCustomGBForceKernel::~ReferenceCalcCustomGBForceKernel() {
    disposeRealArray(particleParamArray, numParticles);
    if (neighborList != NULL)
        delete neighborList;
}

void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
    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.");
        }
    }
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147

    // 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();
    particleParamArray = allocateRealArray(numParticles, numPerParticleParameters);
    for (int i = 0; i < numParticles; ++i) {
        vector<double> parameters;
        force.getParticleParameters(i, parameters);
        for (int j = 0; j < numPerParticleParameters; j++)
            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
    }
    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());
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();

    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
1148
1149
1150
1151
1152
1153
1154
1155
    for (int i = 0; i < force.getNumFunctions(); i++) {
        string name;
        vector<double> values;
        double min, max;
        bool interpolating;
        force.getFunctionParameters(i, name, values, min, max, interpolating);
        functions[name] = new ReferenceTabulatedFunction(min, max, values, interpolating);
    }
1156
1157
1158

    // Parse the expressions for computed values.

1159
    valueDerivExpressions.resize(force.getNumComputedValues());
1160
    valueGradientExpressions.resize(force.getNumComputedValues());
1161
1162
1163
1164
1165
1166
1167
1168
    for (int i = 0; i < force.getNumComputedValues(); i++) {
        string name, expression;
        CustomGBForce::ComputationType type;
        force.getComputedValueParameters(i, name, expression, type);
        Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
        valueExpressions.push_back(ex.createProgram());
        valueTypes.push_back(type);
        valueNames.push_back(name);
1169
1170
1171
        if (i == 0)
            valueDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
        else {
1172
1173
1174
            valueGradientExpressions[i].push_back(ex.differentiate("x").optimize().createProgram());
            valueGradientExpressions[i].push_back(ex.differentiate("y").optimize().createProgram());
            valueGradientExpressions[i].push_back(ex.differentiate("z").optimize().createProgram());
1175
1176
1177
            for (int j = 0; j < i; j++)
                valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
        }
1178
1179
    }

1180
    // Parse the expressions for energy terms.
1181
1182

    energyDerivExpressions.resize(force.getNumEnergyTerms());
1183
    energyGradientExpressions.resize(force.getNumEnergyTerms());
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
    for (int i = 0; i < force.getNumEnergyTerms(); i++) {
        string expression;
        CustomGBForce::ComputationType type;
        force.getEnergyTermParameters(i, expression, type);
        Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
        energyExpressions.push_back(ex.createProgram());
        energyTypes.push_back(type);
        if (type != CustomGBForce::SingleParticle)
            energyDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
        for (int j = 0; j < force.getNumComputedValues(); j++) {
1194
            if (type == CustomGBForce::SingleParticle) {
1195
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
1196
1197
1198
1199
                energyGradientExpressions[i].push_back(ex.differentiate("x").optimize().createProgram());
                energyGradientExpressions[i].push_back(ex.differentiate("y").optimize().createProgram());
                energyGradientExpressions[i].push_back(ex.differentiate("z").optimize().createProgram());
            }
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
            else {
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").optimize().createProgram());
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").optimize().createProgram());
            }
        }
    }

    // Delete the custom functions.

    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
        delete iter->second;
}

void ReferenceCalcCustomGBForceKernel::executeForces(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
1216
1217
    ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
        energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
1218
    bool periodic = (nonbondedMethod == CutoffPeriodic);
1219
1220
    if (periodic)
        ixn.setPeriodic(extractBoxSize(context));
1221
    if (nonbondedMethod != NoCutoff) {
1222
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
1223
1224
1225
1226
1227
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
1228
    ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, 0);
1229
1230
1231
1232
1233
1234
}

double ReferenceCalcCustomGBForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
    RealOpenMM energy = 0;
1235
1236
    ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
        energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
1237
    bool periodic = (nonbondedMethod == CutoffPeriodic);
1238
1239
    if (periodic)
        ixn.setPeriodic(extractBoxSize(context));
1240
    if (nonbondedMethod != NoCutoff) {
1241
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
1242
1243
1244
1245
1246
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
1247
    ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, &energy);
1248
1249
1250
1251
    disposeRealArray(forceData, numParticles);
    return energy;
}

1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
ReferenceCalcCustomExternalForceKernel::~ReferenceCalcCustomExternalForceKernel() {
    disposeRealArray(particleParamArray, numParticles);
}

void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, const CustomExternalForce& force) {
    numParticles = force.getNumParticles();
    int numParameters = force.getNumPerParticleParameters();

    // Build the arrays.

    particles.resize(numParticles);
    particleParamArray = allocateRealArray(numParticles, numParameters);
    vector<double> params;
    for (int i = 0; i < numParticles; ++i) {
        force.getParticleParameters(i, particles[i], params);
        for (int j = 0; j < numParameters; j++)
            particleParamArray[i][j] = (RealOpenMM) params[j];
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
    energyExpression = expression.createProgram();
    forceExpressionX = expression.differentiate("x").optimize().createProgram();
    forceExpressionY = expression.differentiate("y").optimize().createProgram();
    forceExpressionZ = expression.differentiate("z").optimize().createProgram();
    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));
}

void ReferenceCalcCustomExternalForceKernel::executeForces(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceCustomExternalIxn force(energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ, parameterNames, globalParameters);
    for (int i = 0; i < numParticles; ++i)
        force.calculateForce(particles[i], posData, particleParamArray[i], forceData, 0);
}

double ReferenceCalcCustomExternalForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
    RealOpenMM energy = 0;
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceCustomExternalIxn force(energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ, parameterNames, globalParameters);
    for (int i = 0; i < numParticles; ++i)
        force.calculateForce(particles[i], posData, particleParamArray[i], forceData, &energy);
    disposeRealArray(forceData, context.getSystem().getNumParticles());
    return energy;
}

1309
1310
1311
1312
ReferenceCalcCustomHbondForceKernel::~ReferenceCalcCustomHbondForceKernel() {
    disposeRealArray(donorParamArray, numDonors);
    disposeRealArray(acceptorParamArray, numAcceptors);
    disposeIntArray(exclusionArray, numDonors);
1313
1314
    if (ixn != NULL)
        delete ixn;
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
}

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.

1333
    vector<vector<int> > donorParticles(numDonors);
1334
1335
1336
1337
    int numDonorParameters = force.getNumPerDonorParameters();
    donorParamArray = allocateRealArray(numDonors, numDonorParameters);
    for (int i = 0; i < numDonors; ++i) {
        vector<double> parameters;
1338
1339
1340
1341
1342
        int d1, d2, d3;
        force.getDonorParameters(i, d1, d2, d3, parameters);
        donorParticles[i].push_back(d1);
        donorParticles[i].push_back(d2);
        donorParticles[i].push_back(d3);
1343
1344
1345
        for (int j = 0; j < numDonorParameters; j++)
            donorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
    }
1346
    vector<vector<int> > acceptorParticles(numAcceptors);
1347
1348
1349
1350
    int numAcceptorParameters = force.getNumPerAcceptorParameters();
    acceptorParamArray = allocateRealArray(numAcceptors, numAcceptorParameters);
    for (int i = 0; i < numAcceptors; ++i) {
        vector<double> parameters;
1351
1352
1353
1354
1355
        int a1, a2, a3;
        force.getAcceptorParameters(i, a1, a2, a3, parameters);
        acceptorParticles[i].push_back(a1);
        acceptorParticles[i].push_back(a2);
        acceptorParticles[i].push_back(a3);
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
        for (int j = 0; j < numAcceptorParameters; j++)
            acceptorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
    }
    exclusionArray = new int*[numDonors];
    for (int i = 0; i < numDonors; ++i) {
        exclusionArray[i] = new int[exclusions[i].size()+1];
        exclusionArray[i][0] = exclusions[i].size();
        int index = 0;
        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
            exclusionArray[i][++index] = *iter;
    }
1367
    NonbondedMethod nonbondedMethod = CalcCustomHbondForceKernel::NonbondedMethod(force.getNonbondedMethod());
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();

    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
    for (int i = 0; i < force.getNumFunctions(); i++) {
        string name;
        vector<double> values;
        double min, max;
        bool interpolating;
        force.getFunctionParameters(i, name, values, min, max, interpolating);
        functions[name] = new ReferenceTabulatedFunction(min, max, values, interpolating);
    }

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

1384
1385
1386
    map<string, vector<int> > distances;
    map<string, vector<int> > angles;
    map<string, vector<int> > dihedrals;
1387
    Lepton::ParsedExpression energyExpression = CustomHbondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
1388
1389
    vector<string> donorParameterNames;
    vector<string> acceptorParameterNames;
1390
1391
1392
1393
1394
1395
    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));
1396
    ixn = new ReferenceCustomHbondIxn(donorParticles, acceptorParticles, energyExpression, donorParameterNames, acceptorParameterNames, distances, angles, dihedrals);
1397
    isPeriodic = (nonbondedMethod == CutoffPeriodic);
1398
1399
    if (nonbondedMethod != NoCutoff)
        ixn->setUseCutoff(nonbondedCutoff);
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409

    // Delete the custom functions.

    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
        delete iter->second;
}

void ReferenceCalcCustomHbondForceKernel::executeForces(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
1410
1411
    if (isPeriodic)
        ixn->setPeriodic(extractBoxSize(context));
1412
1413
1414
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
1415
    ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, 0);
1416
1417
1418
1419
1420
}

double ReferenceCalcCustomHbondForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
1421
1422
    if (isPeriodic)
        ixn->setPeriodic(extractBoxSize(context));
1423
1424
1425
1426
    RealOpenMM energy = 0;
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
1427
    ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, &energy);
1428
1429
1430
1431
    disposeRealArray(forceData, numParticles);
    return energy;
}

1432
1433
1434
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
    if (dynamics)
        delete dynamics;
1435
1436
    if (constraints)
        delete constraints;
1437
1438
1439
1440
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
1441
1442
    if (constraintDistances)
        delete[] constraintDistances;
1443
1444
}

1445
void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
1446
1447
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
1448
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1449
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1450
1451
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
1452
    constraintDistances = new RealOpenMM[numConstraints];
1453
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
1454
        int particle1, particle2;
1455
        double distance;
Peter Eastman's avatar
Peter Eastman committed
1456
1457
1458
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
1459
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
1460
    }
1461
1462
}

1463
void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
1464
    double stepSize = integrator.getStepSize();
1465
1466
1467
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1468
1469
1470
1471
1472
    if (dynamics == 0 || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
1473
            delete constraints;
1474
        }
Peter Eastman's avatar
Peter Eastman committed
1475
        dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize) );
1476
1477
1478
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
1479
        dynamics->setReferenceConstraintAlgorithm(constraints);
1480
1481
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
1482
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
1483
    data.time += stepSize;
1484
    data.stepCount++;
1485
}
1486

1487
1488
1489
ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
1490
1491
    if (constraints)
        delete constraints;
1492
1493
1494
1495
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
1496
1497
    if (constraintDistances)
        delete[] constraintDistances;
1498
}
1499

1500
void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
1501
1502
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
1503
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1504
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1505
1506
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
1507
    constraintDistances = new RealOpenMM[numConstraints];
1508
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
1509
        int particle1, particle2;
1510
        double distance;
Peter Eastman's avatar
Peter Eastman committed
1511
1512
1513
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
1514
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
1515
    }
1516
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
1517
1518
}

1519
void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
1520
1521
1522
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
1523
1524
1525
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1526
1527
1528
1529
1530
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
1531
            delete constraints;
1532
        }
1533
1534
        RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
        dynamics = new ReferenceStochasticDynamics(
Peter Eastman's avatar
Peter Eastman committed
1535
				context.getSystem().getNumParticles(), 
1536
1537
1538
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(tau), 
				static_cast<RealOpenMM>(temperature) );
1539
1540
1541
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
1542
        dynamics->setReferenceConstraintAlgorithm(constraints);
1543
1544
1545
1546
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
1547
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
1548
    data.time += stepSize;
1549
    data.stepCount++;
1550
1551
}

1552
1553
1554
ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
    if (dynamics)
        delete dynamics;
1555
1556
    if (constraints)
        delete constraints;
1557
1558
1559
1560
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
1561
1562
    if (constraintDistances)
        delete[] constraintDistances;
1563
1564
}

1565
void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
1566
1567
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
1568
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1569
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1570
1571
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
1572
    constraintDistances = new RealOpenMM[numConstraints];
1573
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
1574
        int particle1, particle2;
1575
        double distance;
Peter Eastman's avatar
Peter Eastman committed
1576
1577
1578
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
1579
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
1580
    }
1581
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
1582
1583
}

1584
void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
1585
1586
1587
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
1588
1589
1590
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1591
1592
1593
1594
1595
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
1596
            delete constraints;
1597
        }
1598
        dynamics = new ReferenceBrownianDynamics(
Peter Eastman's avatar
Peter Eastman committed
1599
				context.getSystem().getNumParticles(), 
1600
1601
1602
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(friction), 
				static_cast<RealOpenMM>(temperature) );
1603
1604
1605
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
1606
        dynamics->setReferenceConstraintAlgorithm(constraints);
1607
1608
1609
1610
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
1611
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
1612
    data.time += stepSize;
1613
    data.stepCount++;
1614
1615
}

1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
ReferenceIntegrateVariableLangevinStepKernel::~ReferenceIntegrateVariableLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
    if (constraints)
        delete constraints;
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
    if (constraintDistances)
        delete[] constraintDistances;
}

void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& system, const VariableLangevinIntegrator& integrator) {
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
    for (int i = 0; i < numParticles; ++i)
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
    constraintDistances = new RealOpenMM[numConstraints];
    for (int i = 0; i < numConstraints; ++i) {
        int particle1, particle2;
        double distance;
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
    }
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
}

1648
void ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
1649
1650
1651
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double errorTol = integrator.getErrorTolerance();
1652
1653
1654
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || errorTol != prevErrorTol) {
        // Recreate the computation objects with the new parameters.

        if (dynamics) {
            delete dynamics;
            delete constraints;
        }
        RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
        dynamics = new ReferenceVariableStochasticDynamics(context.getSystem().getNumParticles(), (RealOpenMM) tau, (RealOpenMM) temperature, (RealOpenMM) errorTol);
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
        dynamics->setReferenceConstraintAlgorithm(constraints);
        prevTemp = temperature;
        prevFriction = friction;
        prevErrorTol = errorTol;
    }
    RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
    data.time += dynamics->getDeltaT();
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
}

1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
ReferenceIntegrateVariableVerletStepKernel::~ReferenceIntegrateVariableVerletStepKernel() {
    if (dynamics)
        delete dynamics;
    if (constraints)
        delete constraints;
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
    if (constraintDistances)
        delete[] constraintDistances;
}

void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system, const VariableVerletIntegrator& integrator) {
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
    for (int i = 0; i < numParticles; ++i)
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
    constraintDistances = new RealOpenMM[numConstraints];
    for (int i = 0; i < numConstraints; ++i) {
        int particle1, particle2;
        double distance;
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
    }
}

1711
void ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
1712
    double errorTol = integrator.getErrorTolerance();
1713
1714
1715
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1716
    if (dynamics == 0 || errorTol != prevErrorTol) {
1717
1718
1719
1720
1721
1722
        // Recreate the computation objects with the new parameters.

        if (dynamics) {
            delete dynamics;
            delete constraints;
        }
1723
        dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), (RealOpenMM) errorTol);
1724
1725
1726
1727
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
        dynamics->setReferenceConstraintAlgorithm(constraints);
1728
        prevErrorTol = errorTol;
1729
    }
1730
1731
    RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
1732
    data.time += dynamics->getDeltaT();
1733
1734
1735
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
1736
1737
}

1738
1739
1740
1741
1742
1743
1744
ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
    if (thermostat)
        delete thermostat;
    if (masses)
        delete[] masses;
}

1745
void ReferenceApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
Peter Eastman's avatar
Peter Eastman committed
1746
1747
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
1748
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1749
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1750
    this->thermostat = new ReferenceAndersenThermostat();
1751
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
1752
1753
}

1754
void ReferenceApplyAndersenThermostatKernel::execute(ContextImpl& context) {
1755
    RealOpenMM** velData = extractVelocities(context);
1756
    thermostat->applyThermostat(
1757
			context.getSystem().getNumParticles(),
1758
1759
			velData, 
			masses, 
1760
1761
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::Temperature())), 
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::CollisionFrequency())), 
1762
			static_cast<RealOpenMM>(context.getIntegrator().getStepSize()) );
1763
1764
}

1765
1766
1767
1768
1769
1770
1771
1772
1773
ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel() {
    if (barostat)
        delete barostat;
}

void ReferenceApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& barostat) {
}

void ReferenceApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
1774
1775
    if (barostat == NULL)
        barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), context.getMolecules());
1776
    RealOpenMM** posData = extractPositions(context);
1777
    RealOpenMM* boxSize = extractBoxSize(context);
1778
1779
1780
1781
1782
1783
1784
1785
    barostat->applyBarostat(posData, boxSize, scale);
}

void ReferenceApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    barostat->restorePositions(posData);
}

1786
void ReferenceCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
1787
1788
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
1789
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1790
        masses[i] = system.getParticleMass(i);
1791
1792
}

1793
double ReferenceCalcKineticEnergyKernel::execute(ContextImpl& context) {
1794
    RealOpenMM** velData = extractVelocities(context);
1795
    double energy = 0.0;
1796
    for (size_t i = 0; i < masses.size(); ++i)
1797
1798
        energy += masses[i]*(velData[i][0]*velData[i][0]+velData[i][1]*velData[i][1]+velData[i][2]*velData[i][2]);
    return 0.5*energy;
1799
}
1800

1801
1802
void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
    frequency = force.getFrequency();
Peter Eastman's avatar
Peter Eastman committed
1803
    masses.resize(system.getNumParticles());
1804
    for (size_t i = 0; i < masses.size(); ++i)
Peter Eastman's avatar
Peter Eastman committed
1805
        masses[i] = system.getParticleMass(i);
1806
1807
}

1808
void ReferenceRemoveCMMotionKernel::execute(ContextImpl& context) {
1809
    if (data.stepCount%frequency != 0)
1810
        return;
1811
    RealOpenMM** velData = extractVelocities(context);
1812
1813
1814
1815
    
    // Calculate the center of mass momentum.
    
    RealOpenMM momentum[] = {0.0, 0.0, 0.0};
1816
    RealOpenMM mass = 0.0;
1817
1818
1819
1820
    for (size_t i = 0; i < masses.size(); ++i) {
        momentum[0] += static_cast<RealOpenMM>( masses[i]*velData[i][0] );
        momentum[1] += static_cast<RealOpenMM>( masses[i]*velData[i][1] );
        momentum[2] += static_cast<RealOpenMM>( masses[i]*velData[i][2] );
1821
        mass += static_cast<RealOpenMM>( masses[i] );
1822
1823
    }
    
Peter Eastman's avatar
Peter Eastman committed
1824
    // Adjust the particle velocities.
1825
    
1826
1827
1828
    momentum[0] /= mass;
    momentum[1] /= mass;
    momentum[2] /= mass;
1829
    for (size_t i = 0; i < masses.size(); ++i) {
1830
1831
1832
        velData[i][0] -= momentum[0];
        velData[i][1] -= momentum[1];
        velData[i][2] -= momentum[2];
1833
1834
    }
}