ReferenceKernels.cpp 43.8 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/ReferenceVariableVerletDynamics.h"
37
#include "SimTKReference/ReferenceAndersenThermostat.h"
38
39
#include "SimTKReference/ReferenceAngleBondIxn.h"
#include "SimTKReference/ReferenceBondForce.h"
40
#include "SimTKReference/ReferenceBrownianDynamics.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
#include "SimTKReference/ReferenceCCMAAlgorithm.h"
48
#include "SimTKReference/ReferenceVerletDynamics.h"
49
50
51
52
#include "openmm/CMMotionRemover.h"
#include "openmm/System.h"
#include "openmm/internal/OpenMMContextImpl.h"
#include "openmm/Integrator.h"
53
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
54
#include <cmath>
55
#include <limits>
56
57
58
59

using namespace OpenMM;
using namespace std;

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

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

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

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

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

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

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

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

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

132
133
134
135
136
137
138
139
140
141
142
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;
}

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

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

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

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
217
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
218
219
220
221
222
223
224
    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
225
    disposeRealArray(forceData, context.getSystem().getNumParticles());
226
227
228
229
230
231
232
233
234
235
236
237
238
239
    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
240
        int particle1, particle2, particle3, particle4, periodicity;
241
        double phase, k;
Peter Eastman's avatar
Peter Eastman committed
242
243
244
245
246
        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;
247
248
249
        torsionParamArray[i][0] = (RealOpenMM) k;
        torsionParamArray[i][1] = (RealOpenMM) phase;
        torsionParamArray[i][2] = (RealOpenMM) periodicity;
250
    }
251
252
253
254
255
256
257
258
259
260
261
262
}

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
263
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
264
265
266
267
268
269
270
    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
271
    disposeRealArray(forceData, context.getSystem().getNumParticles());
272
273
274
275
276
277
278
279
280
281
282
283
284
285
    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
286
        int particle1, particle2, particle3, particle4;
287
        double c0, c1, c2, c3, c4, c5;
Peter Eastman's avatar
Peter Eastman committed
288
289
290
291
292
        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;
293
294
295
296
297
298
        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;
299
    }
300
301
302
303
304
305
306
307
308
309
310
311
}

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
312
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
313
314
315
316
317
318
319
    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
320
    disposeRealArray(forceData, context.getSystem().getNumParticles());
321
322
323
324
325
    delete[] energyArray;
    return energy;
}

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

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

    // Identify which exceptions are 1-4 interactions.

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

414
void ReferenceCalcNonbondedForceKernel::executeForces(OpenMMContextImpl& context) {
415
416
    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();
417
    ReferenceLJCoulombIxn clj;
418
    bool periodic = (nonbondedMethod == CutoffPeriodic);
419
    bool ewald  = (nonbondedMethod == Ewald);
420
    bool pme  = (nonbondedMethod == PME);
421
    if (nonbondedMethod != NoCutoff) {
422
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
423
        clj.setUseCutoff(nonbondedCutoff, *neighborList, 78.3f);
424
    }
425
    if (periodic||ewald||pme)
426
        clj.setPeriodic(periodicBoxSize);
427
428
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
429
430
    if (pme)
        clj.setUsePME(ewaldAlpha);
431
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
432
    ReferenceBondForce refBondForce;
433
    ReferenceLJCoulomb14 nonbonded14;
434
    if (nonbondedMethod != NoCutoff)
435
        nonbonded14.setUseCutoff(nonbondedCutoff, 78.3f);
436
437
438
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
}

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

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

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

    // 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
502
                obcParameters->setUseCutoff(static_cast<RealOpenMM>(nonbonded->getCutoffDistance()));
503
504
505
506
507
508
509
510
511
512
513
514
            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;
        }
    }
515
516
    obc = new CpuObc(obcParameters);
    obc->setIncludeAceApproximation(true);
517
518
}

519
void ReferenceCalcGBSAOBCForceKernel::executeForces(OpenMMContextImpl& context) {
520
521
    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();
522
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
523
524
}

525
double ReferenceCalcGBSAOBCForceKernel::executeEnergy(OpenMMContextImpl& context) {
526
    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
527
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
528
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
Peter Eastman's avatar
Peter Eastman committed
529
    disposeRealArray(forceData, context.getSystem().getNumParticles());
530
    return obc->getEnergy();
531
532
}

Mark Friedrichs's avatar
Mark Friedrichs committed
533
534
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
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);
}

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

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

637
638
639
640
641
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
642
643
644
645
646
    if (dynamics == 0 || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
647
            delete constraints;
648
        }
Peter Eastman's avatar
Peter Eastman committed
649
        dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize) );
650
651
652
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
653
        dynamics->setReferenceConstraintAlgorithm(constraints);
654
655
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
656
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
657
    data.time += stepSize;
658
    data.stepCount++;
659
}
660

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

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

693
694
695
696
697
698
699
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
700
701
702
703
704
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
705
            delete constraints;
706
        }
707
708
        RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
        dynamics = new ReferenceStochasticDynamics(
Peter Eastman's avatar
Peter Eastman committed
709
				context.getSystem().getNumParticles(), 
710
711
712
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(tau), 
				static_cast<RealOpenMM>(temperature) );
713
714
715
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
716
        dynamics->setReferenceConstraintAlgorithm(constraints);
717
718
719
720
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
721
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
722
    data.time += stepSize;
723
    data.stepCount++;
724
725
}

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

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

758
759
760
761
762
763
764
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
765
766
767
768
769
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
770
            delete constraints;
771
        }
772
        dynamics = new ReferenceBrownianDynamics(
Peter Eastman's avatar
Peter Eastman committed
773
				context.getSystem().getNumParticles(), 
774
775
776
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(friction), 
				static_cast<RealOpenMM>(temperature) );
777
778
779
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
780
        dynamics->setReferenceConstraintAlgorithm(constraints);
781
782
783
784
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
785
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
786
    data.time += stepSize;
787
    data.stepCount++;
788
789
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
}

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

821
void ReferenceIntegrateVariableVerletStepKernel::execute(OpenMMContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
822
    double errorTol = integrator.getErrorTolerance();
823
824
825
    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
826
    if (dynamics == 0 || errorTol != prevErrorTol) {
827
828
829
830
831
832
        // Recreate the computation objects with the new parameters.

        if (dynamics) {
            delete dynamics;
            delete constraints;
        }
833
        dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), (RealOpenMM) errorTol);
834
835
836
837
        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);
838
        prevErrorTol = errorTol;
839
    }
840
841
    RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
842
    data.time += dynamics->getDeltaT();
843
844
845
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
846
847
}

848
849
850
851
852
853
854
ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
    if (thermostat)
        delete thermostat;
    if (masses)
        delete[] masses;
}

855
void ReferenceApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
Peter Eastman's avatar
Peter Eastman committed
856
857
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
858
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
859
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
860
    this->thermostat = new ReferenceAndersenThermostat();
861
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
862
863
}

864
865
void ReferenceApplyAndersenThermostatKernel::execute(OpenMMContextImpl& context) {
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
866
    thermostat->applyThermostat(
867
			context.getVelocities().getSize(), 
868
869
			velData, 
			masses, 
870
871
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::Temperature())), 
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::CollisionFrequency())), 
872
			static_cast<RealOpenMM>(context.getIntegrator().getStepSize()) );
873
874
}

875
void ReferenceCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
876
877
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
878
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
879
        masses[i] = system.getParticleMass(i);
880
881
}

882
883
double ReferenceCalcKineticEnergyKernel::execute(OpenMMContextImpl& context) {
    RealOpenMM** velData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData()); // Reference code needs to be made const correct
884
    double energy = 0.0;
885
    for (size_t i = 0; i < masses.size(); ++i)
886
887
        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;
888
}
889

890
891
void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
    frequency = force.getFrequency();
Peter Eastman's avatar
Peter Eastman committed
892
    masses.resize(system.getNumParticles());
893
    for (size_t i = 0; i < masses.size(); ++i)
Peter Eastman's avatar
Peter Eastman committed
894
        masses[i] = system.getParticleMass(i);
895
896
}

897
void ReferenceRemoveCMMotionKernel::execute(OpenMMContextImpl& context) {
898
    if (data.stepCount%frequency != 0)
899
900
        return;
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
901
902
903
904
    
    // Calculate the center of mass momentum.
    
    RealOpenMM momentum[] = {0.0, 0.0, 0.0};
905
906
907
908
    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] );
909
910
    }
    
Peter Eastman's avatar
Peter Eastman committed
911
    // Adjust the particle velocities.
912
    
913
914
915
916
917
918
919
    momentum[0] /= static_cast<RealOpenMM>( masses.size() );
    momentum[1] /= static_cast<RealOpenMM>( masses.size() );
    momentum[2] /= static_cast<RealOpenMM>( masses.size() );
    for (size_t i = 0; i < masses.size(); ++i) {
        velData[i][0] -= static_cast<RealOpenMM>( momentum[0]/masses[i] );
        velData[i][1] -= static_cast<RealOpenMM>( momentum[1]/masses[i] );
        velData[i][2] -= static_cast<RealOpenMM>( momentum[2]/masses[i] );
920
921
    }
}