ReferenceKernels.cpp 54.6 KB
Newer Older
1
2
3
4
5
6
7
8
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
9
 * Portions copyright (c) 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/ReferenceCustomNonbondedIxn.h"
41
42
43
44
45
#include "SimTKReference/ReferenceHarmonicBondIxn.h"
#include "SimTKReference/ReferenceLJCoulomb14.h"
#include "SimTKReference/ReferenceLJCoulombIxn.h"
#include "SimTKReference/ReferenceProperDihedralBond.h"
#include "SimTKReference/ReferenceRbDihedralBond.h"
46
#include "SimTKReference/ReferenceStochasticDynamics.h"
47
48
#include "SimTKReference/ReferenceVariableStochasticDynamics.h"
#include "SimTKReference/ReferenceVariableVerletDynamics.h"
49
#include "SimTKReference/ReferenceVerletDynamics.h"
50
51
#include "openmm/CMMotionRemover.h"
#include "openmm/System.h"
52
#include "openmm/internal/ContextImpl.h"
53
#include "openmm/internal/NonbondedForceImpl.h"
54
#include "openmm/Integrator.h"
55
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
56
#include "lepton/CustomFunction.h"
57
58
#include "lepton/Parser.h"
#include "lepton/ParsedExpression.h"
59
#include <cmath>
60
#include <limits>
61
62
63
64

using namespace OpenMM;
using namespace std;

65
static int** allocateIntArray(int length, int width) {
66
67
68
69
70
71
    int** array = new int*[length];
    for (int i = 0; i < length; ++i)
        array[i] = new int[width];
    return array;
}

72
static RealOpenMM** allocateRealArray(int length, int width) {
73
74
75
76
77
78
    RealOpenMM** array = new RealOpenMM*[length];
    for (int i = 0; i < length; ++i)
        array[i] = new RealOpenMM[width];
    return array;
}

79
static int** copyToArray(const vector<vector<int> > vec) {
80
81
82
    if (vec.size() == 0)
        return new int*[0];
    int** array = allocateIntArray(vec.size(), vec[0].size());
83
84
    for (size_t i = 0; i < vec.size(); ++i)
        for (size_t j = 0; j < vec[i].size(); ++j)
85
86
87
88
            array[i][j] = vec[i][j];
    return array;
}

89
static RealOpenMM** copyToArray(const vector<vector<double> > vec) {
90
91
92
    if (vec.size() == 0)
        return new RealOpenMM*[0];
    RealOpenMM** array = allocateRealArray(vec.size(), vec[0].size());
93
94
95
    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]);
96
97
98
    return array;
}

99
static void disposeIntArray(int** array, int size) {
100
101
102
103
104
105
106
    if (array) {
        for (int i = 0; i < size; ++i)
            delete[] array[i];
        delete[] array;
    }
}

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

115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
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;
}

130
static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorithm::AngleInfo>& angles) {
131
132
133
134
135
136
137
    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);
138
                angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, (RealOpenMM)angle));
139
140
141
142
143
            }
        }
    }
}

144
void ReferenceCalcForcesAndEnergyKernel::initialize(const System& system) {
145
146
}

147
void ReferenceCalcForcesAndEnergyKernel::beginForceComputation(ContextImpl& context) {
148
149
150
151
152
153
154
    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;
    }
155
156
}

157
158
159
160
161
162
163
164
165
166
void ReferenceCalcForcesAndEnergyKernel::finishForceComputation(ContextImpl& context) {
}

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

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

167
void ReferenceUpdateStateDataKernel::initialize(const System& system) {
168
169
}

170
double ReferenceUpdateStateDataKernel::getTime(const ContextImpl& context) const {
171
172
173
    return data.time;
}

174
void ReferenceUpdateStateDataKernel::setTime(ContextImpl& context, double time) {
175
176
177
    data.time = time;
}

178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
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) {
        posData[i][0] = positions[i][0];
        posData[i][1] = positions[i][1];
        posData[i][2] = positions[i][2];
    }
}

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) {
        velData[i][0] = velocities[i][0];
        velData[i][1] = velocities[i][1];
        velData[i][2] = velocities[i][2];
    }
}

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

222
ReferenceCalcHarmonicBondForceKernel::~ReferenceCalcHarmonicBondForceKernel() {
223
224
225
226
    disposeIntArray(bondIndexArray, numBonds);
    disposeRealArray(bondParamArray, numBonds);
}

227
void ReferenceCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
228
229
230
231
    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
232
        int particle1, particle2;
233
        double length, k;
Peter Eastman's avatar
Peter Eastman committed
234
235
236
        force.getBondParameters(i, particle1, particle2, length, k);
        bondIndexArray[i][0] = particle1;
        bondIndexArray[i][1] = particle2;
237
238
        bondParamArray[i][0] = (RealOpenMM) length;
        bondParamArray[i][1] = (RealOpenMM) k;
239
    }
240
241
}

242
void ReferenceCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
243
244
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
245
246
247
248
249
    ReferenceBondForce refBondForce;
    ReferenceHarmonicBondIxn harmonicBond;
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, 0, 0, harmonicBond);
}

250
double ReferenceCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context) {
251
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
252
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
253
254
255
256
257
258
259
    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
260
    disposeRealArray(forceData, context.getSystem().getNumParticles());
261
262
263
264
265
266
267
268
269
270
271
272
273
    delete[] energyArray;
    return energy;
}

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);
274
    for (int i = 0; i < force.getNumAngles(); ++i) {
Peter Eastman's avatar
Peter Eastman committed
275
        int particle1, particle2, particle3;
276
        double angle, k;
Peter Eastman's avatar
Peter Eastman committed
277
278
279
280
        force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
        angleIndexArray[i][0] = particle1;
        angleIndexArray[i][1] = particle2;
        angleIndexArray[i][2] = particle3;
281
282
        angleParamArray[i][0] = (RealOpenMM) angle;
        angleParamArray[i][1] = (RealOpenMM) k;
283
    }
284
285
}

286
void ReferenceCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
287
288
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
289
290
291
292
293
    ReferenceBondForce refBondForce;
    ReferenceAngleBondIxn angleBond;
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, angleBond);
}

294
double ReferenceCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context) {
295
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
296
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
297
298
299
300
301
302
303
    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
304
    disposeRealArray(forceData, context.getSystem().getNumParticles());
305
306
307
308
309
310
311
312
313
314
315
316
317
318
    delete[] energyArray;
    return energy;
}

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
319
        int particle1, particle2, particle3, particle4, periodicity;
320
        double phase, k;
Peter Eastman's avatar
Peter Eastman committed
321
322
323
324
325
        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;
326
327
328
        torsionParamArray[i][0] = (RealOpenMM) k;
        torsionParamArray[i][1] = (RealOpenMM) phase;
        torsionParamArray[i][2] = (RealOpenMM) periodicity;
329
    }
330
331
}

332
void ReferenceCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
333
334
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
335
336
337
338
339
    ReferenceBondForce refBondForce;
    ReferenceProperDihedralBond periodicTorsionBond;
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, periodicTorsionBond);
}

340
double ReferenceCalcPeriodicTorsionForceKernel::executeEnergy(ContextImpl& context) {
341
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
342
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
343
344
345
346
347
348
349
    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
350
    disposeRealArray(forceData, context.getSystem().getNumParticles());
351
352
353
354
355
356
357
358
359
360
361
362
363
364
    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
365
        int particle1, particle2, particle3, particle4;
366
        double c0, c1, c2, c3, c4, c5;
Peter Eastman's avatar
Peter Eastman committed
367
368
369
370
371
        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;
372
373
374
375
376
377
        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;
378
    }
379
380
}

381
void ReferenceCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
382
383
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
384
385
386
387
388
    ReferenceBondForce refBondForce;
    ReferenceRbDihedralBond rbTorsionBond;
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, rbTorsionBond);
}

389
double ReferenceCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
390
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
391
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
392
393
394
395
396
397
398
    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
399
    disposeRealArray(forceData, context.getSystem().getNumParticles());
400
401
402
403
404
    delete[] energyArray;
    return energy;
}

ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
Peter Eastman's avatar
Peter Eastman committed
405
406
    disposeRealArray(particleParamArray, numParticles);
    disposeIntArray(exclusionArray, numParticles);
407
408
409
410
411
412
    disposeIntArray(bonded14IndexArray, num14);
    disposeRealArray(bonded14ParamArray, num14);
    if (neighborList != NULL)
        delete neighborList;
}

413
414
415
416
void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {

    // Identify which exceptions are 1-4 interactions.

Peter Eastman's avatar
Peter Eastman committed
417
    numParticles = force.getNumParticles();
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
    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();
433
434
    bonded14IndexArray = allocateIntArray(num14, 2);
    bonded14ParamArray = allocateRealArray(num14, 3);
Peter Eastman's avatar
Peter Eastman committed
435
    particleParamArray = allocateRealArray(numParticles, 3);
436
    RealOpenMM sqrtEps = static_cast<RealOpenMM>( std::sqrt(138.935485) );
Peter Eastman's avatar
Peter Eastman committed
437
    for (int i = 0; i < numParticles; ++i) {
438
        double charge, radius, depth;
Peter Eastman's avatar
Peter Eastman committed
439
440
441
442
        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));
        particleParamArray[i][2] = static_cast<RealOpenMM>(charge*sqrtEps);
443
    }
444
    this->exclusions = exclusions;
Peter Eastman's avatar
Peter Eastman committed
445
446
    exclusionArray = new int*[numParticles];
    for (int i = 0; i < numParticles; ++i) {
447
448
449
450
451
452
453
        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
454
        int particle1, particle2;
455
        double charge, radius, depth;
456
        force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth);
Peter Eastman's avatar
Peter Eastman committed
457
458
        bonded14IndexArray[i][0] = particle1;
        bonded14IndexArray[i][1] = particle2;
459
460
461
        bonded14ParamArray[i][0] = static_cast<RealOpenMM>(radius);
        bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
        bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge*sqrtEps*sqrtEps);
462
    }
463
    nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
464
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
465
    Vec3 boxVectors[3];
466
    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
467
468
469
    periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
    periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
    periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
470
471
472
473
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();
474
475
476
477
478
479
480
481
482
    if (nonbondedMethod == Ewald) {
        double alpha;
        NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmax[0], kmax[1], kmax[2]);
        ewaldAlpha = alpha;
    }
    else if (nonbondedMethod == PME) {
        double alpha;
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2]);
        ewaldAlpha = alpha;
483
    }
484
    rfDielectric = (RealOpenMM)force.getReactionFieldDielectric();
485
486
}

487
void ReferenceCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
488
489
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
490
    ReferenceLJCoulombIxn clj;
491
    bool periodic = (nonbondedMethod == CutoffPeriodic);
492
    bool ewald  = (nonbondedMethod == Ewald);
493
    bool pme  = (nonbondedMethod == PME);
494
    if (nonbondedMethod != NoCutoff) {
495
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
496
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
497
    }
498
    if (periodic||ewald||pme)
499
        clj.setPeriodic(periodicBoxSize);
500
501
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
502
    if (pme)
503
        clj.setUsePME(ewaldAlpha, gridSize);
504
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
505
    ReferenceBondForce refBondForce;
506
    ReferenceLJCoulomb14 nonbonded14;
Peter Eastman's avatar
Peter Eastman committed
507
    if (nonbondedMethod == CutoffNonPeriodic || nonbondedMethod == CutoffPeriodic)
508
        nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
509
510
511
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
}

512
double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
513
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
514
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
515
516
    RealOpenMM energy = 0;
    ReferenceLJCoulombIxn clj;
517
    bool periodic = (nonbondedMethod == CutoffPeriodic);
518
    bool ewald  = (nonbondedMethod == Ewald);
519
    bool pme  = (nonbondedMethod == PME);
520
    if (nonbondedMethod != NoCutoff) {
521
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
522
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
523
    }
524
    if (periodic || ewald || pme)
525
        clj.setPeriodic(periodicBoxSize);
526
527
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
528
    if (pme)
529
        clj.setUsePME(ewaldAlpha, gridSize);
530
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
531
    ReferenceBondForce refBondForce;
532
    ReferenceLJCoulomb14 nonbonded14;
Peter Eastman's avatar
Peter Eastman committed
533
    if (nonbondedMethod == CutoffNonPeriodic || nonbondedMethod == CutoffPeriodic)
534
        nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
Peter Eastman's avatar
Peter Eastman committed
535
536
    RealOpenMM* energyArray = new RealOpenMM[num14];
    for (int i = 0; i < num14; ++i)
537
538
        energyArray[i] = 0;
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, energyArray, 0, &energy, nonbonded14);
Peter Eastman's avatar
Peter Eastman committed
539
    disposeRealArray(forceData, numParticles);
540
541
542
543
    delete[] energyArray;
    return energy;
}

544
545
546
547
548
549
550
551
552
553
554
555
556
class ReferenceCalcCustomNonbondedForceKernel::TabulatedFunction : public Lepton::CustomFunction {
public:
    TabulatedFunction(double min, double max, const vector<double>& values, bool interpolating) :
            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();
557
558
        double scale = (length-1)/(max-min);
        int index = std::floor((x-min)*scale);
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
        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;
        }
576
        x = (x-min)*scale-index;
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
    }
    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);
593
        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.
594
595
596
597
598
599
600
601
602
    }
    CustomFunction* clone() const {
        return new TabulatedFunction(min, max, values, interpolating);
    }
    double min, max;
    vector<double> values;
    bool interpolating;
};

603
604
605
ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
    disposeRealArray(particleParamArray, numParticles);
    disposeIntArray(exclusionArray, numParticles);
606
607
    disposeIntArray(exceptionIndexArray, numExceptions);
    disposeRealArray(exceptionParamArray, numExceptions);
608
609
610
611
612
613
    if (neighborList != NULL)
        delete neighborList;
}

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

614
    // Identify which exceptions are actual interactions.
615
616
617

    numParticles = force.getNumParticles();
    exclusions.resize(numParticles);
618
    vector<int> exceptions;
619
620
621
622
623
624
625
    vector<double> parameters;
    for (int i = 0; i < force.getNumExceptions(); i++) {
        int particle1, particle2;
        force.getExceptionParameters(i, particle1, particle2, parameters);
        exclusions[particle1].insert(particle2);
        exclusions[particle2].insert(particle1);
        if (parameters.size() > 0)
626
            exceptions.push_back(i);
627
628
629
630
    }

    // Build the arrays.

631
    numExceptions = exceptions.size();
632
    int numParameters = force.getNumParameters();
633
634
    exceptionIndexArray = allocateIntArray(numExceptions, 2);
    exceptionParamArray = allocateRealArray(numExceptions, numParameters);
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
    particleParamArray = allocateRealArray(numParticles, numParameters);
    for (int i = 0; i < numParticles; ++i) {
        force.getParticleParameters(i, parameters);
        for (int j = 0; j < numParameters; j++)
            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
    }
    this->exclusions = exclusions;
    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;
    }
650
    for (int i = 0; i < numExceptions; ++i) {
651
        int particle1, particle2;
652
653
654
        force.getExceptionParameters(exceptions[i], particle1, particle2, parameters);
        exceptionIndexArray[i][0] = particle1;
        exceptionIndexArray[i][1] = particle2;
655
        for (int j = 0; j < numParameters; j++)
656
            exceptionParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
657
658
659
660
    }
    nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
    Vec3 boxVectors[3];
661
    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
662
663
664
665
666
667
668
669
    periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
    periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
    periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();

670
671
672
673
674
675
676
677
678
679
680
681
    // 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 TabulatedFunction(min, max, values, interpolating);
    }

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

684
    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
685
686
    energyExpression = expression.createProgram();
    forceExpression = expression.differentiate("r").optimize().createProgram();
687
688
    for (int i = 0; i < numParameters; i++) {
        parameterNames.push_back(force.getParameterName(i));
689
        combiningRules.push_back(Lepton::Parser::parse(force.getParameterCombiningRule(i), functions).optimize().createProgram());
690
691
692
    }
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
693
694
695
696
697

    // Delete the custom functions.

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

void ReferenceCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
701
702
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
703
    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, combiningRules);
704
705
706
707
708
709
710
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    if (periodic)
        ixn.setPeriodic(periodicBoxSize);
711
    map<string, double> globalParameters;
712
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
713
714
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, 0);
715
    ixn.calculateExceptionIxn(numExceptions, exceptionIndexArray, posData, exceptionParamArray, globalParameters, forceData, 0, 0);
716
717
718
}

double ReferenceCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
719
    RealOpenMM** posData = extractPositions(context);
720
721
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
    RealOpenMM energy = 0;
722
    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, combiningRules);
723
724
725
726
727
728
729
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    if (periodic)
        ixn.setPeriodic(periodicBoxSize);
730
    map<string, double> globalParameters;
731
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
732
733
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, &energy);
734
    ixn.calculateExceptionIxn(numExceptions, exceptionIndexArray, posData, exceptionParamArray, globalParameters, forceData, 0, &energy);
735
736
737
738
    disposeRealArray(forceData, numParticles);
    return energy;
}

739
ReferenceCalcGBSAOBCForceKernel::~ReferenceCalcGBSAOBCForceKernel() {
740
    if (obc) {
741
        // delete obc->getObcParameters();
742
743
744
745
        delete obc;
    }
}

746
void ReferenceCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
Peter Eastman's avatar
Peter Eastman committed
747
748
749
750
751
    int numParticles = system.getNumParticles();
    charges.resize(numParticles);
    vector<RealOpenMM> atomicRadii(numParticles);
    vector<RealOpenMM> scaleFactors(numParticles);
    for (int i = 0; i < numParticles; ++i) {
752
        double charge, radius, scalingFactor;
Peter Eastman's avatar
Peter Eastman committed
753
        force.getParticleParameters(i, charge, radius, scalingFactor);
754
755
756
        charges[i] = static_cast<RealOpenMM>(charge);
        atomicRadii[i] = static_cast<RealOpenMM>(radius);
        scaleFactors[i] = static_cast<RealOpenMM>(scalingFactor);
757
    }
Peter Eastman's avatar
Peter Eastman committed
758
    ObcParameters* obcParameters  = new ObcParameters(numParticles, ObcParameters::ObcTypeII);
759
    obcParameters->setAtomicRadii(atomicRadii);
760
    obcParameters->setScaledRadiusFactors(scaleFactors);
761
762
    obcParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
    obcParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
763
764
765
766
767
768
769
770
771
772
    if (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff)
        obcParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
    if (force.getNonbondedMethod() == GBSAOBCForce::CutoffPeriodic) {
        Vec3 boxVectors[3];
        system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
        RealOpenMM periodicBoxSize[3];
        periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
        periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
        periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
        obcParameters->setPeriodic(periodicBoxSize);
773
    }
774
775
    obc = new CpuObc(obcParameters);
    obc->setIncludeAceApproximation(true);
776
777
}

778
void ReferenceCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
779
780
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
781
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
782
783
}

784
double ReferenceCalcGBSAOBCForceKernel::executeEnergy(ContextImpl& context) {
785
    RealOpenMM** posData = extractPositions(context);
Peter Eastman's avatar
Peter Eastman committed
786
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
787
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
Peter Eastman's avatar
Peter Eastman committed
788
    disposeRealArray(forceData, context.getSystem().getNumParticles());
789
    return obc->getEnergy();
790
791
}

Mark Friedrichs's avatar
Mark Friedrichs committed
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
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);
816
817
818
819
820
821
822
823
824
825
826
827
    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()));
    if (force.getNonbondedMethod() == GBVIForce::CutoffPeriodic) {
        Vec3 boxVectors[3];
        system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
        RealOpenMM periodicBoxSize[3];
        periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
        periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
        periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
        gBVIParameters->setPeriodic(periodicBoxSize);
Mark Friedrichs's avatar
Mark Friedrichs committed
828
829
830
831
    }
    gbvi = new CpuGBVI(gBVIParameters);
}

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

834
835
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** forceData = extractForces(context);
Mark Friedrichs's avatar
Mark Friedrichs committed
836
837
838
839
840
841
    RealOpenMM* bornRadii  = new RealOpenMM[context.getSystem().getNumParticles()];
    gbvi->computeBornRadii(posData, bornRadii, NULL ); 
    gbvi->computeBornForces(bornRadii, posData, &charges[0], forceData);
    delete[] bornRadii;
}

842
double ReferenceCalcGBVIForceKernel::executeEnergy(ContextImpl& context) {
843
    RealOpenMM** posData = extractPositions(context);
Mark Friedrichs's avatar
Mark Friedrichs committed
844
845
846
847
848
849
850
    RealOpenMM* bornRadii = new RealOpenMM[context.getSystem().getNumParticles()];
    gbvi->computeBornRadii(posData, bornRadii, NULL ); 
    RealOpenMM energy     = gbvi->computeBornEnergy(bornRadii ,posData, &charges[0]);
    delete[] bornRadii;
    return static_cast<double>(energy);
}

851
852
853
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
    if (dynamics)
        delete dynamics;
854
855
    if (constraints)
        delete constraints;
856
857
858
859
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
860
861
    if (constraintDistances)
        delete[] constraintDistances;
862
863
}

864
void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
865
866
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
867
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
868
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
869
870
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
871
    constraintDistances = new RealOpenMM[numConstraints];
872
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
873
        int particle1, particle2;
874
        double distance;
Peter Eastman's avatar
Peter Eastman committed
875
876
877
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
878
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
879
    }
880
881
}

882
void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
883
    double stepSize = integrator.getStepSize();
884
885
886
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
887
888
889
890
891
    if (dynamics == 0 || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
892
            delete constraints;
893
        }
Peter Eastman's avatar
Peter Eastman committed
894
        dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize) );
895
896
897
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
898
        dynamics->setReferenceConstraintAlgorithm(constraints);
899
900
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
901
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
902
    data.time += stepSize;
903
    data.stepCount++;
904
}
905

906
907
908
ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
909
910
    if (constraints)
        delete constraints;
911
912
913
914
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
915
916
    if (constraintDistances)
        delete[] constraintDistances;
917
}
918

919
void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
920
921
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
922
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
923
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
924
925
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
926
    constraintDistances = new RealOpenMM[numConstraints];
927
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
928
        int particle1, particle2;
929
        double distance;
Peter Eastman's avatar
Peter Eastman committed
930
931
932
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
933
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
934
    }
935
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
936
937
}

938
void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
939
940
941
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
942
943
944
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
945
946
947
948
949
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
950
            delete constraints;
951
        }
952
953
        RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
        dynamics = new ReferenceStochasticDynamics(
Peter Eastman's avatar
Peter Eastman committed
954
				context.getSystem().getNumParticles(), 
955
956
957
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(tau), 
				static_cast<RealOpenMM>(temperature) );
958
959
960
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
961
        dynamics->setReferenceConstraintAlgorithm(constraints);
962
963
964
965
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
966
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
967
    data.time += stepSize;
968
    data.stepCount++;
969
970
}

971
972
973
ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
    if (dynamics)
        delete dynamics;
974
975
    if (constraints)
        delete constraints;
976
977
978
979
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
980
981
    if (constraintDistances)
        delete[] constraintDistances;
982
983
}

984
void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
985
986
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
987
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
988
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
989
990
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
991
    constraintDistances = new RealOpenMM[numConstraints];
992
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
993
        int particle1, particle2;
994
        double distance;
Peter Eastman's avatar
Peter Eastman committed
995
996
997
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
998
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
999
    }
1000
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
1001
1002
}

1003
void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
1004
1005
1006
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
1007
1008
1009
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1010
1011
1012
1013
1014
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
1015
            delete constraints;
1016
        }
1017
        dynamics = new ReferenceBrownianDynamics(
Peter Eastman's avatar
Peter Eastman committed
1018
				context.getSystem().getNumParticles(), 
1019
1020
1021
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(friction), 
				static_cast<RealOpenMM>(temperature) );
1022
1023
1024
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
1025
        dynamics->setReferenceConstraintAlgorithm(constraints);
1026
1027
1028
1029
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
1030
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
1031
    data.time += stepSize;
1032
    data.stepCount++;
1033
1034
}

1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
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());
}

1067
void ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
1068
1069
1070
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double errorTol = integrator.getErrorTolerance();
1071
1072
1073
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
    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++;
}

1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
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);
    }
}

1130
void ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
1131
    double errorTol = integrator.getErrorTolerance();
1132
1133
1134
    RealOpenMM** posData = extractPositions(context);
    RealOpenMM** velData = extractVelocities(context);
    RealOpenMM** forceData = extractForces(context);
1135
    if (dynamics == 0 || errorTol != prevErrorTol) {
1136
1137
1138
1139
1140
1141
        // Recreate the computation objects with the new parameters.

        if (dynamics) {
            delete dynamics;
            delete constraints;
        }
1142
        dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), (RealOpenMM) errorTol);
1143
1144
1145
1146
        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);
1147
        prevErrorTol = errorTol;
1148
    }
1149
1150
    RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
1151
    data.time += dynamics->getDeltaT();
1152
1153
1154
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
1155
1156
}

1157
1158
1159
1160
1161
1162
1163
ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
    if (thermostat)
        delete thermostat;
    if (masses)
        delete[] masses;
}

1164
void ReferenceApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
Peter Eastman's avatar
Peter Eastman committed
1165
1166
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
1167
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1168
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1169
    this->thermostat = new ReferenceAndersenThermostat();
1170
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
1171
1172
}

1173
void ReferenceApplyAndersenThermostatKernel::execute(ContextImpl& context) {
1174
    RealOpenMM** velData = extractVelocities(context);
1175
    thermostat->applyThermostat(
1176
			context.getSystem().getNumParticles(),
1177
1178
			velData, 
			masses, 
1179
1180
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::Temperature())), 
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::CollisionFrequency())), 
1181
			static_cast<RealOpenMM>(context.getIntegrator().getStepSize()) );
1182
1183
}

1184
void ReferenceCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
1185
1186
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
1187
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1188
        masses[i] = system.getParticleMass(i);
1189
1190
}

1191
double ReferenceCalcKineticEnergyKernel::execute(ContextImpl& context) {
1192
    RealOpenMM** velData = extractVelocities(context);
1193
    double energy = 0.0;
1194
    for (size_t i = 0; i < masses.size(); ++i)
1195
1196
        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;
1197
}
1198

1199
1200
void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
    frequency = force.getFrequency();
Peter Eastman's avatar
Peter Eastman committed
1201
    masses.resize(system.getNumParticles());
1202
    for (size_t i = 0; i < masses.size(); ++i)
Peter Eastman's avatar
Peter Eastman committed
1203
        masses[i] = system.getParticleMass(i);
1204
1205
}

1206
void ReferenceRemoveCMMotionKernel::execute(ContextImpl& context) {
1207
    if (data.stepCount%frequency != 0)
1208
        return;
1209
    RealOpenMM** velData = extractVelocities(context);
1210
1211
1212
1213
    
    // Calculate the center of mass momentum.
    
    RealOpenMM momentum[] = {0.0, 0.0, 0.0};
1214
    RealOpenMM mass = 0.0;
1215
1216
1217
1218
    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] );
1219
        mass += static_cast<RealOpenMM>( masses[i] );
1220
1221
    }
    
Peter Eastman's avatar
Peter Eastman committed
1222
    // Adjust the particle velocities.
1223
    
1224
1225
1226
    momentum[0] /= mass;
    momentum[1] /= mass;
    momentum[2] /= mass;
1227
    for (size_t i = 0; i < masses.size(); ++i) {
1228
1229
1230
        velData[i][0] -= momentum[0];
        velData[i][1] -= momentum[1];
        velData[i][2] -= momentum[2];
1231
1232
    }
}