ReferenceKernels.cpp 80.1 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/ReferenceCustomAngleIxn.h"
41
#include "SimTKReference/ReferenceCustomBondIxn.h"
42
#include "SimTKReference/ReferenceCustomExternalIxn.h"
43
#include "SimTKReference/ReferenceCustomGBIxn.h"
44
#include "SimTKReference/ReferenceCustomHbondIxn.h"
45
#include "SimTKReference/ReferenceCustomNonbondedIxn.h"
46
#include "SimTKReference/ReferenceCustomTorsionIxn.h"
47
48
49
#include "SimTKReference/ReferenceHarmonicBondIxn.h"
#include "SimTKReference/ReferenceLJCoulomb14.h"
#include "SimTKReference/ReferenceLJCoulombIxn.h"
50
#include "SimTKReference/ReferenceMonteCarloBarostat.h"
51
52
#include "SimTKReference/ReferenceProperDihedralBond.h"
#include "SimTKReference/ReferenceRbDihedralBond.h"
53
#include "SimTKReference/ReferenceStochasticDynamics.h"
54
55
#include "SimTKReference/ReferenceVariableStochasticDynamics.h"
#include "SimTKReference/ReferenceVariableVerletDynamics.h"
56
#include "SimTKReference/ReferenceVerletDynamics.h"
57
#include "openmm/CMMotionRemover.h"
58
#include "openmm/Context.h"
59
#include "openmm/System.h"
60
#include "openmm/internal/ContextImpl.h"
61
#include "openmm/internal/CustomHbondForceImpl.h"
62
#include "openmm/internal/NonbondedForceImpl.h"
63
#include "openmm/Integrator.h"
64
#include "openmm/OpenMMException.h"
65
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
66
#include "lepton/CustomFunction.h"
67
68
#include "lepton/Parser.h"
#include "lepton/ParsedExpression.h"
69
#include <cmath>
70
#include <limits>
71
72
73
74

using namespace OpenMM;
using namespace std;

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

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

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

99
static RealOpenMM** copyToArray(const vector<vector<double> > vec) {
100
101
102
    if (vec.size() == 0)
        return new RealOpenMM*[0];
    RealOpenMM** array = allocateRealArray(vec.size(), vec[0].size());
103
104
105
    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]);
106
107
108
    return array;
}

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

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

125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
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;
}

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

145
static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorithm::AngleInfo>& angles) {
146
147
148
149
150
151
152
    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);
153
                angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, (RealOpenMM)angle));
154
155
156
157
158
            }
        }
    }
}

159
void ReferenceCalcForcesAndEnergyKernel::initialize(const System& system) {
160
161
}

162
void ReferenceCalcForcesAndEnergyKernel::beginForceComputation(ContextImpl& context) {
163
164
165
166
167
168
169
    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;
    }
170
171
}

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

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

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

182
void ReferenceUpdateStateDataKernel::initialize(const System& system) {
183
184
}

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

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

193
194
195
196
197
198
199
200
201
202
203
204
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) {
205
206
207
        posData[i][0] = (RealOpenMM) positions[i][0];
        posData[i][1] = (RealOpenMM) positions[i][1];
        posData[i][2] = (RealOpenMM) positions[i][2];
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
    }
}

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) {
223
224
225
        velData[i][0] = (RealOpenMM) velocities[i][0];
        velData[i][1] = (RealOpenMM) velocities[i][1];
        velData[i][2] = (RealOpenMM) velocities[i][2];
226
227
228
229
230
231
232
233
234
235
236
    }
}

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]);
}

237
238
239
240
241
242
243
244
245
246
247
248
249
250
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];
}

251
ReferenceCalcHarmonicBondForceKernel::~ReferenceCalcHarmonicBondForceKernel() {
252
253
254
255
    disposeIntArray(bondIndexArray, numBonds);
    disposeRealArray(bondParamArray, numBonds);
}

256
void ReferenceCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
257
258
259
260
    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
261
        int particle1, particle2;
262
        double length, k;
Peter Eastman's avatar
Peter Eastman committed
263
264
265
        force.getBondParameters(i, particle1, particle2, length, k);
        bondIndexArray[i][0] = particle1;
        bondIndexArray[i][1] = particle2;
266
267
        bondParamArray[i][0] = (RealOpenMM) length;
        bondParamArray[i][1] = (RealOpenMM) k;
268
    }
269
270
}

271
void ReferenceCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
272
273
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
274
275
276
277
278
    ReferenceBondForce refBondForce;
    ReferenceHarmonicBondIxn harmonicBond;
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, 0, 0, harmonicBond);
}

279
double ReferenceCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
280
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
281
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
282
283
284
285
286
287
288
    RealOpenMM* energyArray = new RealOpenMM[numBonds];
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceHarmonicBondIxn harmonicBond;
    for (int i = 0; i < numBonds; ++i)
        energyArray[i] = 0;
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, energyArray, 0, &energy, harmonicBond);
Peter Eastman's avatar
Peter Eastman committed
289
    disposeRealArray(forceData, context.getSystem().getNumParticles());
290
291
292
293
    delete[] energyArray;
    return energy;
}

294
295
296
297
298
299
300
301
302
303
304
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.

305
    bondIndexArray = allocateIntArray(numBonds, 2);
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
    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);
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, 0, 0, harmonicBond);
}

double ReferenceCalcCustomBondForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
    RealOpenMM* energyArray = new RealOpenMM[numBonds];
    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);
    for (int i = 0; i < numBonds; ++i)
        energyArray[i] = 0;
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, energyArray, 0, &energy, harmonicBond);
    disposeRealArray(forceData, context.getSystem().getNumParticles());
    delete[] energyArray;
    return energy;
}

357
358
359
360
361
362
363
364
365
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);
366
    for (int i = 0; i < force.getNumAngles(); ++i) {
Peter Eastman's avatar
Peter Eastman committed
367
        int particle1, particle2, particle3;
368
        double angle, k;
Peter Eastman's avatar
Peter Eastman committed
369
370
371
372
        force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
        angleIndexArray[i][0] = particle1;
        angleIndexArray[i][1] = particle2;
        angleIndexArray[i][2] = particle3;
373
374
        angleParamArray[i][0] = (RealOpenMM) angle;
        angleParamArray[i][1] = (RealOpenMM) k;
375
    }
376
377
}

378
void ReferenceCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
379
380
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
381
382
383
384
385
    ReferenceBondForce refBondForce;
    ReferenceAngleBondIxn angleBond;
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, angleBond);
}

386
double ReferenceCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
387
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
388
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
389
390
391
392
393
394
395
    RealOpenMM* energyArray = new RealOpenMM[numAngles];
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceAngleBondIxn angleBond;
    for (int i = 0; i < numAngles; ++i)
        energyArray[i] = 0;
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, energyArray, 0, &energy, angleBond);
Peter Eastman's avatar
Peter Eastman committed
396
    disposeRealArray(forceData, context.getSystem().getNumParticles());
397
398
399
400
    delete[] energyArray;
    return energy;
}

401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
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;
443
444
    ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters);
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, customAngle);
445
446
447
448
449
450
451
452
453
454
455
}

double ReferenceCalcCustomAngleForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
    RealOpenMM* energyArray = new RealOpenMM[numAngles];
    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;
456
    ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters);
457
458
    for (int i = 0; i < numAngles; ++i)
        energyArray[i] = 0;
459
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, energyArray, 0, &energy, customAngle);
460
461
462
463
464
    disposeRealArray(forceData, context.getSystem().getNumParticles());
    delete[] energyArray;
    return energy;
}

465
466
467
468
469
470
471
472
473
474
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
475
        int particle1, particle2, particle3, particle4, periodicity;
476
        double phase, k;
Peter Eastman's avatar
Peter Eastman committed
477
478
479
480
481
        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;
482
483
484
        torsionParamArray[i][0] = (RealOpenMM) k;
        torsionParamArray[i][1] = (RealOpenMM) phase;
        torsionParamArray[i][2] = (RealOpenMM) periodicity;
485
    }
486
487
}

488
void ReferenceCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
489
490
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
491
492
493
494
495
    ReferenceBondForce refBondForce;
    ReferenceProperDihedralBond periodicTorsionBond;
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, periodicTorsionBond);
}

496
double ReferenceCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
497
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
498
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
499
500
501
502
503
504
505
    RealOpenMM* energyArray = new RealOpenMM[numTorsions];
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceProperDihedralBond periodicTorsionBond;
    for (int i = 0; i < numTorsions; ++i)
        energyArray[i] = 0;
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, energyArray, 0, &energy, periodicTorsionBond);
Peter Eastman's avatar
Peter Eastman committed
506
    disposeRealArray(forceData, context.getSystem().getNumParticles());
507
508
509
510
511
512
513
514
515
516
517
518
519
520
    delete[] energyArray;
    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
521
        int particle1, particle2, particle3, particle4;
522
        double c0, c1, c2, c3, c4, c5;
Peter Eastman's avatar
Peter Eastman committed
523
524
525
526
527
        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;
528
529
530
531
532
533
        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;
534
    }
535
536
}

537
void ReferenceCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
538
539
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
540
541
542
543
544
    ReferenceBondForce refBondForce;
    ReferenceRbDihedralBond rbTorsionBond;
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, rbTorsionBond);
}

545
double ReferenceCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
546
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
547
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
548
549
550
551
552
553
554
    RealOpenMM* energyArray = new RealOpenMM[numTorsions];
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceRbDihedralBond rbTorsionBond;
    for (int i = 0; i < numTorsions; ++i)
        energyArray[i] = 0;
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, energyArray, 0, &energy, rbTorsionBond);
Peter Eastman's avatar
Peter Eastman committed
555
    disposeRealArray(forceData, context.getSystem().getNumParticles());
556
557
558
559
    delete[] energyArray;
    return energy;
}

560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
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
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);
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, customTorsion);
}

double ReferenceCalcCustomTorsionForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
    RealOpenMM* energyArray = new RealOpenMM[numTorsions];
    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);
    for (int i = 0; i < numTorsions; ++i)
        energyArray[i] = 0;
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, energyArray, 0, &energy, customTorsion);
    disposeRealArray(forceData, context.getSystem().getNumParticles());
    delete[] energyArray;
    return energy;
}

625
ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
Peter Eastman's avatar
Peter Eastman committed
626
627
    disposeRealArray(particleParamArray, numParticles);
    disposeIntArray(exclusionArray, numParticles);
628
629
630
631
632
633
    disposeIntArray(bonded14IndexArray, num14);
    disposeRealArray(bonded14ParamArray, num14);
    if (neighborList != NULL)
        delete neighborList;
}

634
635
636
637
void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {

    // Identify which exceptions are 1-4 interactions.

Peter Eastman's avatar
Peter Eastman committed
638
    numParticles = force.getNumParticles();
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
    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();
654
655
    bonded14IndexArray = allocateIntArray(num14, 2);
    bonded14ParamArray = allocateRealArray(num14, 3);
Peter Eastman's avatar
Peter Eastman committed
656
657
    particleParamArray = allocateRealArray(numParticles, 3);
    for (int i = 0; i < numParticles; ++i) {
658
        double charge, radius, depth;
Peter Eastman's avatar
Peter Eastman committed
659
660
661
        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));
662
        particleParamArray[i][2] = static_cast<RealOpenMM>(charge);
663
    }
664
    this->exclusions = exclusions;
Peter Eastman's avatar
Peter Eastman committed
665
666
    exclusionArray = new int*[numParticles];
    for (int i = 0; i < numParticles; ++i) {
667
668
669
670
671
672
673
        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
674
        int particle1, particle2;
675
        double charge, radius, depth;
676
        force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth);
Peter Eastman's avatar
Peter Eastman committed
677
678
        bonded14IndexArray[i][0] = particle1;
        bonded14IndexArray[i][1] = particle2;
679
680
        bonded14ParamArray[i][0] = static_cast<RealOpenMM>(radius);
        bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
681
        bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge);
682
    }
683
    nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
684
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
685
686
687
688
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();
689
690
691
    if (nonbondedMethod == Ewald) {
        double alpha;
        NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmax[0], kmax[1], kmax[2]);
692
        ewaldAlpha = (RealOpenMM) alpha;
693
694
695
696
    }
    else if (nonbondedMethod == PME) {
        double alpha;
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2]);
697
        ewaldAlpha = (RealOpenMM) alpha;
698
    }
699
    rfDielectric = (RealOpenMM)force.getReactionFieldDielectric();
700
701
}

702
void ReferenceCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
703
704
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
705
    ReferenceLJCoulombIxn clj;
706
    bool periodic = (nonbondedMethod == CutoffPeriodic);
707
    bool ewald  = (nonbondedMethod == Ewald);
708
    bool pme  = (nonbondedMethod == PME);
709
    if (nonbondedMethod != NoCutoff) {
710
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
711
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
712
    }
713
    if (periodic || ewald || pme)
714
        clj.setPeriodic(extractBoxSize(context));
715
716
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
717
    if (pme)
718
        clj.setUsePME(ewaldAlpha, gridSize);
719
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
720
    ReferenceBondForce refBondForce;
721
722
723
724
    ReferenceLJCoulomb14 nonbonded14;
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
}

725
double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
726
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
727
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
728
729
    RealOpenMM energy = 0;
    ReferenceLJCoulombIxn clj;
730
    bool periodic = (nonbondedMethod == CutoffPeriodic);
731
    bool ewald  = (nonbondedMethod == Ewald);
732
    bool pme  = (nonbondedMethod == PME);
733
    if (nonbondedMethod != NoCutoff) {
734
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
735
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
736
    }
737
    if (periodic || ewald || pme)
738
        clj.setPeriodic(extractBoxSize(context));
739
740
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
741
    if (pme)
742
        clj.setUsePME(ewaldAlpha, gridSize);
743
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
744
    ReferenceBondForce refBondForce;
745
    ReferenceLJCoulomb14 nonbonded14;
Peter Eastman's avatar
Peter Eastman committed
746
747
    RealOpenMM* energyArray = new RealOpenMM[num14];
    for (int i = 0; i < num14; ++i)
748
749
        energyArray[i] = 0;
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, energyArray, 0, &energy, nonbonded14);
Peter Eastman's avatar
Peter Eastman committed
750
    disposeRealArray(forceData, numParticles);
751
752
753
754
    delete[] energyArray;
    return energy;
}

755
class ReferenceTabulatedFunction : public Lepton::CustomFunction {
756
public:
757
    ReferenceTabulatedFunction(double min, double max, const vector<double>& values, bool interpolating) :
758
759
760
761
762
763
764
765
766
767
            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();
768
        double scale = (length-1)/(max-min);
769
        int index = (int) std::floor((x-min)*scale);
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
        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;
        }
787
        x = (x-min)*scale-index;
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
    }
    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);
804
        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.
805
806
    }
    CustomFunction* clone() const {
807
        return new ReferenceTabulatedFunction(min, max, values, interpolating);
808
809
810
811
812
813
    }
    double min, max;
    vector<double> values;
    bool interpolating;
};

814
815
816
817
818
819
820
821
822
ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
    disposeRealArray(particleParamArray, numParticles);
    disposeIntArray(exclusionArray, numParticles);
    if (neighborList != NULL)
        delete neighborList;
}

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

823
    // Record the exclusions.
824
825
826

    numParticles = force.getNumParticles();
    exclusions.resize(numParticles);
827
    for (int i = 0; i < force.getNumExclusions(); i++) {
828
        int particle1, particle2;
829
        force.getExclusionParticles(i, particle1, particle2);
830
831
832
833
834
835
        exclusions[particle1].insert(particle2);
        exclusions[particle2].insert(particle1);
    }

    // Build the arrays.

836
    int numParameters = force.getNumPerParticleParameters();
837
838
    particleParamArray = allocateRealArray(numParticles, numParameters);
    for (int i = 0; i < numParticles; ++i) {
839
        vector<double> parameters;
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
        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();

859
860
861
862
863
864
865
866
867
    // 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);
868
        functions[name] = new ReferenceTabulatedFunction(min, max, values, interpolating);
869
870
    }

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

873
    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
874
875
    energyExpression = expression.createProgram();
    forceExpression = expression.differentiate("r").optimize().createProgram();
876
877
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerParticleParameterName(i));
878
879
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
880
881
882
883
884

    // Delete the custom functions.

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

void ReferenceCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
888
889
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
890
    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames);
891
892
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
893
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
894
895
896
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    if (periodic)
897
        ixn.setPeriodic(extractBoxSize(context));
898
    map<string, double> globalParameters;
899
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
900
901
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, 0);
902
903
904
}

double ReferenceCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
905
    RealOpenMM** posData = extractPositions(context);
906
907
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
    RealOpenMM energy = 0;
908
    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames);
909
910
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
911
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
912
913
914
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    if (periodic)
915
        ixn.setPeriodic(extractBoxSize(context));
916
    map<string, double> globalParameters;
917
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
918
919
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, &energy);
920
921
922
923
    disposeRealArray(forceData, numParticles);
    return energy;
}

924
ReferenceCalcGBSAOBCForceKernel::~ReferenceCalcGBSAOBCForceKernel() {
925
    if (obc) {
926
        // delete obc->getObcParameters();
927
928
929
930
        delete obc;
    }
}

931
void ReferenceCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
Peter Eastman's avatar
Peter Eastman committed
932
933
934
935
936
    int numParticles = system.getNumParticles();
    charges.resize(numParticles);
    vector<RealOpenMM> atomicRadii(numParticles);
    vector<RealOpenMM> scaleFactors(numParticles);
    for (int i = 0; i < numParticles; ++i) {
937
        double charge, radius, scalingFactor;
Peter Eastman's avatar
Peter Eastman committed
938
        force.getParticleParameters(i, charge, radius, scalingFactor);
939
940
941
        charges[i] = static_cast<RealOpenMM>(charge);
        atomicRadii[i] = static_cast<RealOpenMM>(radius);
        scaleFactors[i] = static_cast<RealOpenMM>(scalingFactor);
942
    }
943
    ObcParameters* obcParameters = new ObcParameters(numParticles, ObcParameters::ObcTypeII);
944
    obcParameters->setAtomicRadii(atomicRadii);
945
    obcParameters->setScaledRadiusFactors(scaleFactors);
946
947
    obcParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
    obcParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
948
949
    if (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff)
        obcParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
950
    isPeriodic = (force.getNonbondedMethod() == GBSAOBCForce::CutoffPeriodic);
951
952
    obc = new CpuObc(obcParameters);
    obc->setIncludeAceApproximation(true);
953
954
}

955
void ReferenceCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
956
957
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
958
959
    if (isPeriodic)
        obc->getObcParameters()->setPeriodic(extractBoxSize(context));
960
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
961
962
}

963
double ReferenceCalcGBSAOBCForceKernel::executeEnergy(ContextImpl& context) {
964
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
965
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
966
967
    if (isPeriodic)
        obc->getObcParameters()->setPeriodic(extractBoxSize(context));
968
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
Peter Eastman's avatar
Peter Eastman committed
969
    disposeRealArray(forceData, context.getSystem().getNumParticles());
970
    return obc->getEnergy();
971
972
}

Mark Friedrichs's avatar
Mark Friedrichs committed
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
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);
997
998
999
1000
    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()));
1001
    isPeriodic = (force.getNonbondedMethod() == GBVIForce::CutoffPeriodic);
Mark Friedrichs's avatar
Mark Friedrichs committed
1002
1003
1004
    gbvi = new CpuGBVI(gBVIParameters);
}

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

1007
1008
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
Mark Friedrichs's avatar
Mark Friedrichs committed
1009
    RealOpenMM* bornRadii  = new RealOpenMM[context.getSystem().getNumParticles()];
1010
1011
    if (isPeriodic)
        gbvi->getGBVIParameters()->setPeriodic(extractBoxSize(context));
Mark Friedrichs's avatar
Mark Friedrichs committed
1012
1013
1014
1015
1016
    gbvi->computeBornRadii(posData, bornRadii, NULL ); 
    gbvi->computeBornForces(bornRadii, posData, &charges[0], forceData);
    delete[] bornRadii;
}

1017
double ReferenceCalcGBVIForceKernel::executeEnergy(ContextImpl& context) {
1018
    RealOpenMM** posData = extractPositions(context);
Mark Friedrichs's avatar
Mark Friedrichs committed
1019
    RealOpenMM* bornRadii = new RealOpenMM[context.getSystem().getNumParticles()];
1020
1021
    if (isPeriodic)
        gbvi->getGBVIParameters()->setPeriodic(extractBoxSize(context));
1022
    gbvi->computeBornRadii(posData, bornRadii, NULL );
Mark Friedrichs's avatar
Mark Friedrichs committed
1023
1024
1025
1026
1027
    RealOpenMM energy     = gbvi->computeBornEnergy(bornRadii ,posData, &charges[0]);
    delete[] bornRadii;
    return static_cast<double>(energy);
}

1028
1029
1030
1031
1032
1033
1034
ReferenceCalcCustomGBForceKernel::~ReferenceCalcCustomGBForceKernel() {
    disposeRealArray(particleParamArray, numParticles);
    if (neighborList != NULL)
        delete neighborList;
}

void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
    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.");
        }
    }
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082

    // 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;
1083
1084
1085
1086
1087
1088
1089
1090
    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);
    }
1091
1092
1093

    // Parse the expressions for computed values.

1094
    valueDerivExpressions.resize(force.getNumComputedValues());
1095
    valueGradientExpressions.resize(force.getNumComputedValues());
1096
1097
1098
1099
1100
1101
1102
1103
    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);
1104
1105
1106
        if (i == 0)
            valueDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
        else {
1107
1108
1109
            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());
1110
1111
1112
            for (int j = 0; j < i; j++)
                valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
        }
1113
1114
    }

1115
    // Parse the expressions for energy terms.
1116
1117

    energyDerivExpressions.resize(force.getNumEnergyTerms());
1118
    energyGradientExpressions.resize(force.getNumEnergyTerms());
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
    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++) {
1129
            if (type == CustomGBForce::SingleParticle) {
1130
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
1131
1132
1133
1134
                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());
            }
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
            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);
1151
1152
    ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
        energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
1153
    bool periodic = (nonbondedMethod == CutoffPeriodic);
1154
1155
    if (periodic)
        ixn.setPeriodic(extractBoxSize(context));
1156
    if (nonbondedMethod != NoCutoff) {
1157
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
1158
1159
1160
1161
1162
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
1163
    ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, 0);
1164
1165
1166
1167
1168
1169
}

double ReferenceCalcCustomGBForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
    RealOpenMM energy = 0;
1170
1171
    ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
        energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
1172
    bool periodic = (nonbondedMethod == CutoffPeriodic);
1173
1174
    if (periodic)
        ixn.setPeriodic(extractBoxSize(context));
1175
    if (nonbondedMethod != NoCutoff) {
1176
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? extractBoxSize(context) : NULL, nonbondedCutoff, 0.0);
1177
1178
1179
1180
1181
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
1182
    ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, &energy);
1183
1184
1185
1186
    disposeRealArray(forceData, numParticles);
    return energy;
}

1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
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;
}

1244
1245
1246
1247
ReferenceCalcCustomHbondForceKernel::~ReferenceCalcCustomHbondForceKernel() {
    disposeRealArray(donorParamArray, numDonors);
    disposeRealArray(acceptorParamArray, numAcceptors);
    disposeIntArray(exclusionArray, numDonors);
1248
1249
    if (ixn != NULL)
        delete ixn;
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
}

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.

1268
    vector<vector<int> > donorParticles(numDonors);
1269
1270
1271
1272
    int numDonorParameters = force.getNumPerDonorParameters();
    donorParamArray = allocateRealArray(numDonors, numDonorParameters);
    for (int i = 0; i < numDonors; ++i) {
        vector<double> parameters;
1273
1274
1275
1276
1277
        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);
1278
1279
1280
        for (int j = 0; j < numDonorParameters; j++)
            donorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
    }
1281
    vector<vector<int> > acceptorParticles(numAcceptors);
1282
1283
1284
1285
    int numAcceptorParameters = force.getNumPerAcceptorParameters();
    acceptorParamArray = allocateRealArray(numAcceptors, numAcceptorParameters);
    for (int i = 0; i < numAcceptors; ++i) {
        vector<double> parameters;
1286
1287
1288
1289
1290
        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);
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
        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;
    }
1302
    NonbondedMethod nonbondedMethod = CalcCustomHbondForceKernel::NonbondedMethod(force.getNonbondedMethod());
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
    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);
    }

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

1319
1320
1321
    map<string, vector<int> > distances;
    map<string, vector<int> > angles;
    map<string, vector<int> > dihedrals;
1322
    Lepton::ParsedExpression energyExpression = CustomHbondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
1323
1324
    vector<string> donorParameterNames;
    vector<string> acceptorParameterNames;
1325
1326
1327
1328
1329
1330
    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));
1331
    ixn = new ReferenceCustomHbondIxn(donorParticles, acceptorParticles, energyExpression, donorParameterNames, acceptorParameterNames, distances, angles, dihedrals);
1332
    isPeriodic = (nonbondedMethod == CutoffPeriodic);
1333
1334
    if (nonbondedMethod != NoCutoff)
        ixn->setUseCutoff(nonbondedCutoff);
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344

    // 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);
1345
1346
    if (isPeriodic)
        ixn->setPeriodic(extractBoxSize(context));
1347
1348
1349
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
1350
    ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, 0);
1351
1352
1353
1354
1355
}

double ReferenceCalcCustomHbondForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
1356
1357
    if (isPeriodic)
        ixn->setPeriodic(extractBoxSize(context));
1358
1359
1360
1361
    RealOpenMM energy = 0;
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
1362
    ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, &energy);
1363
1364
1365
1366
    disposeRealArray(forceData, numParticles);
    return energy;
}

1367
1368
1369
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
    if (dynamics)
        delete dynamics;
1370
1371
    if (constraints)
        delete constraints;
1372
1373
1374
1375
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
1376
1377
    if (constraintDistances)
        delete[] constraintDistances;
1378
1379
}

1380
void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
1381
1382
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
1383
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1384
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1385
1386
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
1387
    constraintDistances = new RealOpenMM[numConstraints];
1388
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
1389
        int particle1, particle2;
1390
        double distance;
Peter Eastman's avatar
Peter Eastman committed
1391
1392
1393
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
1394
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
1395
    }
1396
1397
}

1398
void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
1399
    double stepSize = integrator.getStepSize();
1400
1401
1402
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1403
1404
1405
1406
1407
    if (dynamics == 0 || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
1408
            delete constraints;
1409
        }
Peter Eastman's avatar
Peter Eastman committed
1410
        dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize) );
1411
1412
1413
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
1414
        dynamics->setReferenceConstraintAlgorithm(constraints);
1415
1416
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
1417
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
1418
    data.time += stepSize;
1419
    data.stepCount++;
1420
}
1421

1422
1423
1424
ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
1425
1426
    if (constraints)
        delete constraints;
1427
1428
1429
1430
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
1431
1432
    if (constraintDistances)
        delete[] constraintDistances;
1433
}
1434

1435
void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
1436
1437
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
1438
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1439
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1440
1441
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
1442
    constraintDistances = new RealOpenMM[numConstraints];
1443
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
1444
        int particle1, particle2;
1445
        double distance;
Peter Eastman's avatar
Peter Eastman committed
1446
1447
1448
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
1449
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
1450
    }
1451
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
1452
1453
}

1454
void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
1455
1456
1457
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
1458
1459
1460
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1461
1462
1463
1464
1465
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
1466
            delete constraints;
1467
        }
1468
1469
        RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
        dynamics = new ReferenceStochasticDynamics(
Peter Eastman's avatar
Peter Eastman committed
1470
				context.getSystem().getNumParticles(), 
1471
1472
1473
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(tau), 
				static_cast<RealOpenMM>(temperature) );
1474
1475
1476
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
1477
        dynamics->setReferenceConstraintAlgorithm(constraints);
1478
1479
1480
1481
        prevTemp = temperature;
        prevFriction = friction;
        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
ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
    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 ReferenceIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& 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 ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& 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
        dynamics = new ReferenceBrownianDynamics(
Peter Eastman's avatar
Peter Eastman committed
1534
				context.getSystem().getNumParticles(), 
1535
1536
1537
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(friction), 
				static_cast<RealOpenMM>(temperature) );
1538
1539
1540
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
1541
        dynamics->setReferenceConstraintAlgorithm(constraints);
1542
1543
1544
1545
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
1546
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
1547
    data.time += stepSize;
1548
    data.stepCount++;
1549
1550
}

1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
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());
}

1583
void ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
1584
1585
1586
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double errorTol = integrator.getErrorTolerance();
1587
1588
1589
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
    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++;
}

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
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);
    }
}

1646
void ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
1647
    double errorTol = integrator.getErrorTolerance();
1648
1649
1650
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1651
    if (dynamics == 0 || errorTol != prevErrorTol) {
1652
1653
1654
1655
1656
1657
        // Recreate the computation objects with the new parameters.

        if (dynamics) {
            delete dynamics;
            delete constraints;
        }
1658
        dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), (RealOpenMM) errorTol);
1659
1660
1661
1662
        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);
1663
        prevErrorTol = errorTol;
1664
    }
1665
1666
    RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
1667
    data.time += dynamics->getDeltaT();
1668
1669
1670
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
1671
1672
}

1673
1674
1675
1676
1677
1678
1679
ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
    if (thermostat)
        delete thermostat;
    if (masses)
        delete[] masses;
}

1680
void ReferenceApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
Peter Eastman's avatar
Peter Eastman committed
1681
1682
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
1683
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1684
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1685
    this->thermostat = new ReferenceAndersenThermostat();
1686
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
1687
1688
}

1689
void ReferenceApplyAndersenThermostatKernel::execute(ContextImpl& context) {
1690
    RealOpenMM** velData = extractVelocities(context);
1691
    thermostat->applyThermostat(
1692
			context.getSystem().getNumParticles(),
1693
1694
			velData, 
			masses, 
1695
1696
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::Temperature())), 
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::CollisionFrequency())), 
1697
			static_cast<RealOpenMM>(context.getIntegrator().getStepSize()) );
1698
1699
}

1700
1701
1702
1703
1704
1705
1706
1707
1708
ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel() {
    if (barostat)
        delete barostat;
}

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

void ReferenceApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
1709
1710
    if (barostat == NULL)
        barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), context.getMolecules());
1711
    RealOpenMM** posData = extractPositions(context);
1712
    RealOpenMM* boxSize = extractBoxSize(context);
1713
1714
1715
1716
1717
1718
1719
1720
    barostat->applyBarostat(posData, boxSize, scale);
}

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

1721
void ReferenceCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
1722
1723
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
1724
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1725
        masses[i] = system.getParticleMass(i);
1726
1727
}

1728
double ReferenceCalcKineticEnergyKernel::execute(ContextImpl& context) {
1729
    RealOpenMM** velData = extractVelocities(context);
1730
    double energy = 0.0;
1731
    for (size_t i = 0; i < masses.size(); ++i)
1732
1733
        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;
1734
}
1735

1736
1737
void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
    frequency = force.getFrequency();
Peter Eastman's avatar
Peter Eastman committed
1738
    masses.resize(system.getNumParticles());
1739
    for (size_t i = 0; i < masses.size(); ++i)
Peter Eastman's avatar
Peter Eastman committed
1740
        masses[i] = system.getParticleMass(i);
1741
1742
}

1743
void ReferenceRemoveCMMotionKernel::execute(ContextImpl& context) {
1744
    if (data.stepCount%frequency != 0)
1745
        return;
1746
    RealOpenMM** velData = extractVelocities(context);
1747
1748
1749
1750
    
    // Calculate the center of mass momentum.
    
    RealOpenMM momentum[] = {0.0, 0.0, 0.0};
1751
    RealOpenMM mass = 0.0;
1752
1753
1754
1755
    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] );
1756
        mass += static_cast<RealOpenMM>( masses[i] );
1757
1758
    }
    
Peter Eastman's avatar
Peter Eastman committed
1759
    // Adjust the particle velocities.
1760
    
1761
1762
1763
    momentum[0] /= mass;
    momentum[1] /= mass;
    momentum[2] /= mass;
1764
    for (size_t i = 0; i < masses.size(); ++i) {
1765
1766
1767
        velData[i][0] -= momentum[0];
        velData[i][1] -= momentum[1];
        velData[i][2] -= momentum[2];
1768
1769
    }
}