ReferenceKernels.cpp 47.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 "ReferenceFloatStreamImpl.h"
34
#include "gbsa/CpuObc.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
35
#include "gbsa/CpuGBVI.h"
36
#include "SimTKReference/ReferenceAndersenThermostat.h"
37
38
#include "SimTKReference/ReferenceAngleBondIxn.h"
#include "SimTKReference/ReferenceBondForce.h"
39
#include "SimTKReference/ReferenceBrownianDynamics.h"
40
#include "SimTKReference/ReferenceCCMAAlgorithm.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
52
53
#include "openmm/CMMotionRemover.h"
#include "openmm/System.h"
#include "openmm/internal/OpenMMContextImpl.h"
#include "openmm/Integrator.h"
54
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
55
#include <cmath>
56
#include <limits>
57
58
59
60

using namespace OpenMM;
using namespace std;

61
static int** allocateIntArray(int length, int width) {
62
63
64
65
66
67
    int** array = new int*[length];
    for (int i = 0; i < length; ++i)
        array[i] = new int[width];
    return array;
}

68
static RealOpenMM** allocateRealArray(int length, int width) {
69
70
71
72
73
74
    RealOpenMM** array = new RealOpenMM*[length];
    for (int i = 0; i < length; ++i)
        array[i] = new RealOpenMM[width];
    return array;
}

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

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

95
static void disposeIntArray(int** array, int size) {
96
97
98
99
100
101
102
    if (array) {
        for (int i = 0; i < size; ++i)
            delete[] array[i];
        delete[] array;
    }
}

103
static void disposeRealArray(RealOpenMM** array, int size) {
104
105
106
107
108
109
110
    if (array) {
        for (int i = 0; i < size; ++i)
            delete[] array[i];
        delete[] array;
    }
}

111
static void findAnglesForCCMA(const System& system, vector<ReferenceCCMAAlgorithm::AngleInfo>& angles) {
112
113
114
115
116
117
118
    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);
119
                angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, (RealOpenMM)angle));
120
121
122
123
124
            }
        }
    }
}

125
126
127
void ReferenceInitializeForcesKernel::initialize(const System& system) {
}

128
void ReferenceInitializeForcesKernel::execute(OpenMMContextImpl& context) {
129
130
131
132
    double zero[] = {0.0, 0.0, 0.0};
    context.getForces().fillWithValue(zero);
}

133
134
135
136
137
138
139
140
141
142
143
void ReferenceUpdateTimeKernel::initialize(const System& system) {
}

double ReferenceUpdateTimeKernel::getTime(const OpenMMContextImpl& context) const {
    return data.time;
}

void ReferenceUpdateTimeKernel::setTime(OpenMMContextImpl& context, double time) {
    data.time = time;
}

144
ReferenceCalcHarmonicBondForceKernel::~ReferenceCalcHarmonicBondForceKernel() {
145
146
147
148
    disposeIntArray(bondIndexArray, numBonds);
    disposeRealArray(bondParamArray, numBonds);
}

149
void ReferenceCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
150
151
152
153
    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
154
        int particle1, particle2;
155
        double length, k;
Peter Eastman's avatar
Peter Eastman committed
156
157
158
        force.getBondParameters(i, particle1, particle2, length, k);
        bondIndexArray[i][0] = particle1;
        bondIndexArray[i][1] = particle2;
159
160
        bondParamArray[i][0] = (RealOpenMM) length;
        bondParamArray[i][1] = (RealOpenMM) k;
161
    }
162
163
164
165
166
167
168
169
170
171
172
173
}

void ReferenceCalcHarmonicBondForceKernel::executeForces(OpenMMContextImpl& context) {
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
    RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
    ReferenceBondForce refBondForce;
    ReferenceHarmonicBondIxn harmonicBond;
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, 0, 0, harmonicBond);
}

double ReferenceCalcHarmonicBondForceKernel::executeEnergy(OpenMMContextImpl& context) {
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
Peter Eastman's avatar
Peter Eastman committed
174
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
175
176
177
178
179
180
181
    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
182
    disposeRealArray(forceData, context.getSystem().getNumParticles());
183
184
185
186
187
188
189
190
191
192
193
194
195
    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);
196
    for (int i = 0; i < force.getNumAngles(); ++i) {
Peter Eastman's avatar
Peter Eastman committed
197
        int particle1, particle2, particle3;
198
        double angle, k;
Peter Eastman's avatar
Peter Eastman committed
199
200
201
202
        force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
        angleIndexArray[i][0] = particle1;
        angleIndexArray[i][1] = particle2;
        angleIndexArray[i][2] = particle3;
203
204
        angleParamArray[i][0] = (RealOpenMM) angle;
        angleParamArray[i][1] = (RealOpenMM) k;
205
    }
206
207
208
209
210
211
212
213
214
215
216
217
}

void ReferenceCalcHarmonicAngleForceKernel::executeForces(OpenMMContextImpl& context) {
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
    RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
    ReferenceBondForce refBondForce;
    ReferenceAngleBondIxn angleBond;
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, angleBond);
}

double ReferenceCalcHarmonicAngleForceKernel::executeEnergy(OpenMMContextImpl& context) {
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
Peter Eastman's avatar
Peter Eastman committed
218
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
219
220
221
222
223
224
225
    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
226
    disposeRealArray(forceData, context.getSystem().getNumParticles());
227
228
229
230
231
232
233
234
235
236
237
238
239
240
    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
241
        int particle1, particle2, particle3, particle4, periodicity;
242
        double phase, k;
Peter Eastman's avatar
Peter Eastman committed
243
244
245
246
247
        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;
248
249
250
        torsionParamArray[i][0] = (RealOpenMM) k;
        torsionParamArray[i][1] = (RealOpenMM) phase;
        torsionParamArray[i][2] = (RealOpenMM) periodicity;
251
    }
252
253
254
255
256
257
258
259
260
261
262
263
}

void ReferenceCalcPeriodicTorsionForceKernel::executeForces(OpenMMContextImpl& context) {
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
    RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
    ReferenceBondForce refBondForce;
    ReferenceProperDihedralBond periodicTorsionBond;
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, periodicTorsionBond);
}

double ReferenceCalcPeriodicTorsionForceKernel::executeEnergy(OpenMMContextImpl& context) {
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
Peter Eastman's avatar
Peter Eastman committed
264
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
265
266
267
268
269
270
271
    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
272
    disposeRealArray(forceData, context.getSystem().getNumParticles());
273
274
275
276
277
278
279
280
281
282
283
284
285
286
    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
287
        int particle1, particle2, particle3, particle4;
288
        double c0, c1, c2, c3, c4, c5;
Peter Eastman's avatar
Peter Eastman committed
289
290
291
292
293
        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;
294
295
296
297
298
299
        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;
300
    }
301
302
303
304
305
306
307
308
309
310
311
312
}

void ReferenceCalcRBTorsionForceKernel::executeForces(OpenMMContextImpl& context) {
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
    RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
    ReferenceBondForce refBondForce;
    ReferenceRbDihedralBond rbTorsionBond;
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, rbTorsionBond);
}

double ReferenceCalcRBTorsionForceKernel::executeEnergy(OpenMMContextImpl& context) {
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
Peter Eastman's avatar
Peter Eastman committed
313
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
314
315
316
317
318
319
320
    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
321
    disposeRealArray(forceData, context.getSystem().getNumParticles());
322
323
324
325
326
    delete[] energyArray;
    return energy;
}

ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
Peter Eastman's avatar
Peter Eastman committed
327
328
    disposeRealArray(particleParamArray, numParticles);
    disposeIntArray(exclusionArray, numParticles);
329
330
331
332
333
334
    disposeIntArray(bonded14IndexArray, num14);
    disposeRealArray(bonded14ParamArray, num14);
    if (neighborList != NULL)
        delete neighborList;
}

335
336
337
338
void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {

    // Identify which exceptions are 1-4 interactions.

Peter Eastman's avatar
Peter Eastman committed
339
    numParticles = force.getNumParticles();
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
    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();
355
356
    bonded14IndexArray = allocateIntArray(num14, 2);
    bonded14ParamArray = allocateRealArray(num14, 3);
Peter Eastman's avatar
Peter Eastman committed
357
    particleParamArray = allocateRealArray(numParticles, 3);
358
    RealOpenMM sqrtEps = static_cast<RealOpenMM>( std::sqrt(138.935485) );
Peter Eastman's avatar
Peter Eastman committed
359
    for (int i = 0; i < numParticles; ++i) {
360
        double charge, radius, depth;
Peter Eastman's avatar
Peter Eastman committed
361
362
363
364
        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);
365
    }
366
    this->exclusions = exclusions;
Peter Eastman's avatar
Peter Eastman committed
367
368
    exclusionArray = new int*[numParticles];
    for (int i = 0; i < numParticles; ++i) {
369
370
371
372
373
374
375
        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
376
        int particle1, particle2;
377
        double charge, radius, depth;
378
        force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth);
Peter Eastman's avatar
Peter Eastman committed
379
380
        bonded14IndexArray[i][0] = particle1;
        bonded14IndexArray[i][1] = particle2;
381
382
383
        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);
384
    }
385
    nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
386
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
387
388
389
390
391
    Vec3 boxVectors[3];
    force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
    periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
    periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
    periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
392
393
394
395
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();
396
    if (nonbondedMethod == Ewald || nonbondedMethod == PME) {
397
398
399
400
401
402
        RealOpenMM ewaldErrorTol = (RealOpenMM) force.getEwaldErrorTolerance();
        ewaldAlpha = (RealOpenMM) (std::sqrt(-std::log(ewaldErrorTol))/nonbondedCutoff);
        RealOpenMM mx = periodicBoxSize[0]/nonbondedCutoff;
        RealOpenMM my = periodicBoxSize[1]/nonbondedCutoff;
        RealOpenMM mz = periodicBoxSize[2]/nonbondedCutoff;
        RealOpenMM pi = (RealOpenMM) 3.1415926535897932385;
403
404
405
        kmax[0] = (int)std::ceil(-(mx/pi)*std::log(ewaldErrorTol));
        kmax[1] = (int)std::ceil(-(my/pi)*std::log(ewaldErrorTol));
        kmax[2] = (int)std::ceil(-(mz/pi)*std::log(ewaldErrorTol));
406
407
408
409
410
411
        if (kmax[0]%2 == 0)
            kmax[0]++;
        if (kmax[1]%2 == 0)
            kmax[1]++;
        if (kmax[2]%2 == 0)
            kmax[2]++;
412
    }
413
    rfDielectric = force.getReactionFieldDielectric();
414
415
}

416
void ReferenceCalcNonbondedForceKernel::executeForces(OpenMMContextImpl& context) {
417
418
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
    RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
419
    ReferenceLJCoulombIxn clj;
420
    bool periodic = (nonbondedMethod == CutoffPeriodic);
421
    bool ewald  = (nonbondedMethod == Ewald);
422
    bool pme  = (nonbondedMethod == PME);
423
    if (nonbondedMethod != NoCutoff) {
424
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
425
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
426
    }
427
    if (periodic||ewald||pme)
428
        clj.setPeriodic(periodicBoxSize);
429
430
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
431
432
    if (pme)
        clj.setUsePME(ewaldAlpha);
433
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
434
    ReferenceBondForce refBondForce;
435
    ReferenceLJCoulomb14 nonbonded14;
436
    if (nonbondedMethod != NoCutoff)
437
        nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
438
439
440
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
}

441
double ReferenceCalcNonbondedForceKernel::executeEnergy(OpenMMContextImpl& context) {
442
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
Peter Eastman's avatar
Peter Eastman committed
443
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
444
445
    RealOpenMM energy = 0;
    ReferenceLJCoulombIxn clj;
446
    bool periodic = (nonbondedMethod == CutoffPeriodic);
447
    bool ewald  = (nonbondedMethod == Ewald);
448
    bool pme  = (nonbondedMethod == PME);
449
    if (nonbondedMethod != NoCutoff) {
450
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
451
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
452
    }
453
    if (periodic || ewald || pme)
454
        clj.setPeriodic(periodicBoxSize);
455
456
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
457
458
    if (pme)
        clj.setUsePME(ewaldAlpha);
459
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
460
    ReferenceBondForce refBondForce;
461
    ReferenceLJCoulomb14 nonbonded14;
462
    if (nonbondedMethod != NoCutoff)
463
        nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
Peter Eastman's avatar
Peter Eastman committed
464
465
    RealOpenMM* energyArray = new RealOpenMM[num14];
    for (int i = 0; i < num14; ++i)
466
467
        energyArray[i] = 0;
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, energyArray, 0, &energy, nonbonded14);
Peter Eastman's avatar
Peter Eastman committed
468
    disposeRealArray(forceData, numParticles);
469
470
471
472
    delete[] energyArray;
    return energy;
}

473
ReferenceCalcGBSAOBCForceKernel::~ReferenceCalcGBSAOBCForceKernel() {
474
    if (obc) {
475
        // delete obc->getObcParameters();
476
477
478
479
        delete obc;
    }
}

480
void ReferenceCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
Peter Eastman's avatar
Peter Eastman committed
481
482
483
484
485
    int numParticles = system.getNumParticles();
    charges.resize(numParticles);
    vector<RealOpenMM> atomicRadii(numParticles);
    vector<RealOpenMM> scaleFactors(numParticles);
    for (int i = 0; i < numParticles; ++i) {
486
        double charge, radius, scalingFactor;
Peter Eastman's avatar
Peter Eastman committed
487
        force.getParticleParameters(i, charge, radius, scalingFactor);
488
489
490
        charges[i] = static_cast<RealOpenMM>(charge);
        atomicRadii[i] = static_cast<RealOpenMM>(radius);
        scaleFactors[i] = static_cast<RealOpenMM>(scalingFactor);
491
    }
Peter Eastman's avatar
Peter Eastman committed
492
    ObcParameters* obcParameters  = new ObcParameters(numParticles, ObcParameters::ObcTypeII);
493
    obcParameters->setAtomicRadii(atomicRadii);
494
    obcParameters->setScaledRadiusFactors(scaleFactors);
495
496
    obcParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
    obcParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
497
498
499
500
501
502
503

    // If there is a NonbondedForce in this system, use it to initialize cutoffs and periodic boundary conditions.

    for (int i = 0; i < system.getNumForces(); i++) {
        const NonbondedForce* nonbonded = dynamic_cast<const NonbondedForce*>(&system.getForce(i));
        if (nonbonded != NULL) {
            if (nonbonded->getNonbondedMethod() != NonbondedForce::NoCutoff)
Mark Friedrichs's avatar
Mark Friedrichs committed
504
                obcParameters->setUseCutoff(static_cast<RealOpenMM>(nonbonded->getCutoffDistance()));
505
506
507
508
509
510
511
512
513
514
515
516
            if (nonbonded->getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
                Vec3 boxVectors[3];
                nonbonded->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);
            }
            break;
        }
    }
517
518
    obc = new CpuObc(obcParameters);
    obc->setIncludeAceApproximation(true);
519
520
}

521
void ReferenceCalcGBSAOBCForceKernel::executeForces(OpenMMContextImpl& context) {
522
523
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
    RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
524
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
525
526
}

527
double ReferenceCalcGBSAOBCForceKernel::executeEnergy(OpenMMContextImpl& context) {
528
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
Peter Eastman's avatar
Peter Eastman committed
529
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
530
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
Peter Eastman's avatar
Peter Eastman committed
531
    disposeRealArray(forceData, context.getSystem().getNumParticles());
532
    return obc->getEnergy();
533
534
}

Mark Friedrichs's avatar
Mark Friedrichs committed
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
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
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);
    gBVIParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
    gBVIParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );

    // If there is a NonbondedForce in this system, use it to initialize cutoffs and periodic boundary conditions.

    for (int i = 0; i < system.getNumForces(); i++) {
        const NonbondedForce* nonbonded = dynamic_cast<const NonbondedForce*>(&system.getForce(i));
        if (nonbonded != NULL) {
            if (nonbonded->getNonbondedMethod() != NonbondedForce::NoCutoff)
                gBVIParameters->setUseCutoff( static_cast<RealOpenMM>(nonbonded->getCutoffDistance()));
            if (nonbonded->getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
                Vec3 boxVectors[3];
                nonbonded->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);
            }
            break;
        }
    }
    gbvi = new CpuGBVI(gBVIParameters);

}

void ReferenceCalcGBVIForceKernel::executeForces(OpenMMContextImpl& context) {

    RealOpenMM** posData   = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
    RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
    RealOpenMM* bornRadii  = new RealOpenMM[context.getSystem().getNumParticles()];
    gbvi->computeBornRadii(posData, bornRadii, NULL ); 
    gbvi->computeBornForces(bornRadii, posData, &charges[0], forceData);
    delete[] bornRadii;
}

double ReferenceCalcGBVIForceKernel::executeEnergy(OpenMMContextImpl& context) {
    RealOpenMM** posData  = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
    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);
}

608
609
610
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
    if (dynamics)
        delete dynamics;
611
612
    if (constraints)
        delete constraints;
613
614
615
616
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
617
618
    if (constraintDistances)
        delete[] constraintDistances;
619
620
}

621
void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
622
623
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
624
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
625
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
626
627
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
628
    constraintDistances = new RealOpenMM[numConstraints];
629
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
630
        int particle1, particle2;
631
        double distance;
Peter Eastman's avatar
Peter Eastman committed
632
633
634
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
635
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
636
    }
637
638
}

639
640
641
642
643
void ReferenceIntegrateVerletStepKernel::execute(OpenMMContextImpl& context, const VerletIntegrator& integrator) {
    double stepSize = integrator.getStepSize();
    RealOpenMM** posData = ((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData();
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
    RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData()); // Reference code needs to be made const correct
644
645
646
647
648
    if (dynamics == 0 || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
649
            delete constraints;
650
        }
Peter Eastman's avatar
Peter Eastman committed
651
        dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize) );
652
653
654
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
655
        dynamics->setReferenceConstraintAlgorithm(constraints);
656
657
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
658
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
659
    data.time += stepSize;
660
    data.stepCount++;
661
}
662

663
664
665
ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
666
667
    if (constraints)
        delete constraints;
668
669
670
671
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
672
673
    if (constraintDistances)
        delete[] constraintDistances;
674
}
675

676
void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
677
678
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
679
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
680
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
681
682
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
683
    constraintDistances = new RealOpenMM[numConstraints];
684
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
685
        int particle1, particle2;
686
        double distance;
Peter Eastman's avatar
Peter Eastman committed
687
688
689
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
690
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
691
    }
692
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
693
694
}

695
696
697
698
699
700
701
void ReferenceIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, const LangevinIntegrator& integrator) {
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
    RealOpenMM** posData = ((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData();
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
    RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData()); // Reference code needs to be made const correct
702
703
704
705
706
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
707
            delete constraints;
708
        }
709
710
        RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
        dynamics = new ReferenceStochasticDynamics(
Peter Eastman's avatar
Peter Eastman committed
711
				context.getSystem().getNumParticles(), 
712
713
714
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(tau), 
				static_cast<RealOpenMM>(temperature) );
715
716
717
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
718
        dynamics->setReferenceConstraintAlgorithm(constraints);
719
720
721
722
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
723
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
724
    data.time += stepSize;
725
    data.stepCount++;
726
727
}

728
729
730
ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
    if (dynamics)
        delete dynamics;
731
732
    if (constraints)
        delete constraints;
733
734
735
736
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
737
738
    if (constraintDistances)
        delete[] constraintDistances;
739
740
}

741
void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
742
743
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
744
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
745
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
746
747
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
748
    constraintDistances = new RealOpenMM[numConstraints];
749
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
750
        int particle1, particle2;
751
        double distance;
Peter Eastman's avatar
Peter Eastman committed
752
753
754
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
755
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
756
    }
757
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
758
759
}

760
761
762
763
764
765
766
void ReferenceIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, const BrownianIntegrator& integrator) {
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
    RealOpenMM** posData = ((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData();
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
    RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData()); // Reference code needs to be made const correct
767
768
769
770
771
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
772
            delete constraints;
773
        }
774
        dynamics = new ReferenceBrownianDynamics(
Peter Eastman's avatar
Peter Eastman committed
775
				context.getSystem().getNumParticles(), 
776
777
778
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(friction), 
				static_cast<RealOpenMM>(temperature) );
779
780
781
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
782
        dynamics->setReferenceConstraintAlgorithm(constraints);
783
784
785
786
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
787
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
788
    data.time += stepSize;
789
    data.stepCount++;
790
791
}

792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
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());
}

void ReferenceIntegrateVariableLangevinStepKernel::execute(OpenMMContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double errorTol = integrator.getErrorTolerance();
    RealOpenMM** posData = ((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData();
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
    RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData()); // Reference code needs to be made const correct
    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++;
}

856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
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);
    }
}

887
void ReferenceIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
888
    double errorTol = integrator.getErrorTolerance();
889
890
891
    RealOpenMM** posData = ((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData();
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
    RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData()); // Reference code needs to be made const correct
892
    if (dynamics == 0 || errorTol != prevErrorTol) {
893
894
895
896
897
898
        // Recreate the computation objects with the new parameters.

        if (dynamics) {
            delete dynamics;
            delete constraints;
        }
899
        dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), (RealOpenMM) errorTol);
900
901
902
903
        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);
904
        prevErrorTol = errorTol;
905
    }
906
907
    RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
908
    data.time += dynamics->getDeltaT();
909
910
911
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
912
913
}

914
915
916
917
918
919
920
ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
    if (thermostat)
        delete thermostat;
    if (masses)
        delete[] masses;
}

921
void ReferenceApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
Peter Eastman's avatar
Peter Eastman committed
922
923
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
924
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
925
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
926
    this->thermostat = new ReferenceAndersenThermostat();
927
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
928
929
}

930
931
void ReferenceApplyAndersenThermostatKernel::execute(OpenMMContextImpl& context) {
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
932
    thermostat->applyThermostat(
933
			context.getVelocities().getSize(), 
934
935
			velData, 
			masses, 
936
937
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::Temperature())), 
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::CollisionFrequency())), 
938
			static_cast<RealOpenMM>(context.getIntegrator().getStepSize()) );
939
940
}

941
void ReferenceCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
942
943
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
944
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
945
        masses[i] = system.getParticleMass(i);
946
947
}

948
949
double ReferenceCalcKineticEnergyKernel::execute(OpenMMContextImpl& context) {
    RealOpenMM** velData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData()); // Reference code needs to be made const correct
950
    double energy = 0.0;
951
    for (size_t i = 0; i < masses.size(); ++i)
952
953
        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;
954
}
955

956
957
void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
    frequency = force.getFrequency();
Peter Eastman's avatar
Peter Eastman committed
958
    masses.resize(system.getNumParticles());
959
    for (size_t i = 0; i < masses.size(); ++i)
Peter Eastman's avatar
Peter Eastman committed
960
        masses[i] = system.getParticleMass(i);
961
962
}

963
void ReferenceRemoveCMMotionKernel::execute(OpenMMContextImpl& context) {
964
    if (data.stepCount%frequency != 0)
965
966
        return;
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
967
968
969
970
    
    // Calculate the center of mass momentum.
    
    RealOpenMM momentum[] = {0.0, 0.0, 0.0};
971
    RealOpenMM mass = 0.0;
972
973
974
975
    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] );
976
        mass += masses[i];
977
978
    }
    
Peter Eastman's avatar
Peter Eastman committed
979
    // Adjust the particle velocities.
980
    
981
982
983
    momentum[0] /= mass;
    momentum[1] /= mass;
    momentum[2] /= mass;
984
    for (size_t i = 0; i < masses.size(); ++i) {
985
986
987
        velData[i][0] -= momentum[0];
        velData[i][1] -= momentum[1];
        velData[i][2] -= momentum[2];
988
989
    }
}