ReferenceKernels.cpp 24.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
/* -------------------------------------------------------------------------- *
 *                                   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.               *
 *                                                                            *
 * Portions copyright (c) 2008 Stanford University and the Authors.           *
 * 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"
35
#include "SimTKReference/ReferenceAndersenThermostat.h"
36
37
#include "SimTKReference/ReferenceAngleBondIxn.h"
#include "SimTKReference/ReferenceBondForce.h"
38
#include "SimTKReference/ReferenceBrownianDynamics.h"
39
40
41
42
43
#include "SimTKReference/ReferenceHarmonicBondIxn.h"
#include "SimTKReference/ReferenceLJCoulomb14.h"
#include "SimTKReference/ReferenceLJCoulombIxn.h"
#include "SimTKReference/ReferenceProperDihedralBond.h"
#include "SimTKReference/ReferenceRbDihedralBond.h"
44
45
#include "SimTKReference/ReferenceStochasticDynamics.h"
#include "SimTKReference/ReferenceShakeAlgorithm.h"
46
#include "SimTKReference/ReferenceVerletDynamics.h"
47
#include <cmath>
48
#include <limits>
49
50
51
52

using namespace OpenMM;
using namespace std;

53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
int** allocateIntArray(int length, int width) {
    int** array = new int*[length];
    for (int i = 0; i < length; ++i)
        array[i] = new int[width];
    return array;
}

RealOpenMM** allocateRealArray(int length, int width) {
    RealOpenMM** array = new RealOpenMM*[length];
    for (int i = 0; i < length; ++i)
        array[i] = new RealOpenMM[width];
    return array;
}

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

RealOpenMM** copyToArray(const vector<vector<double> > vec) {
    if (vec.size() == 0)
        return new RealOpenMM*[0];
    RealOpenMM** array = allocateRealArray(vec.size(), vec[0].size());
81
82
83
    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]);
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
    return array;
}

void disposeIntArray(int** array, int size) {
    if (array) {
        for (int i = 0; i < size; ++i)
            delete[] array[i];
        delete[] array;
    }
}

void disposeRealArray(RealOpenMM** array, int size) {
    if (array) {
        for (int i = 0; i < size; ++i)
            delete[] array[i];
        delete[] array;
    }
}

103
ReferenceCalcStandardMMForceFieldKernel::~ReferenceCalcStandardMMForceFieldKernel() {
104
105
106
107
108
109
110
111
    disposeIntArray(bondIndexArray, numBonds);
    disposeRealArray(bondParamArray, numBonds);
    disposeIntArray(angleIndexArray, numAngles);
    disposeRealArray(angleParamArray, numAngles);
    disposeIntArray(periodicTorsionIndexArray, numPeriodicTorsions);
    disposeRealArray(periodicTorsionParamArray, numPeriodicTorsions);
    disposeIntArray(rbTorsionIndexArray, numRBTorsions);
    disposeRealArray(rbTorsionParamArray, numRBTorsions);
112
113
114
115
    disposeRealArray(atomParamArray, numAtoms);
    disposeIntArray(exclusionArray, numAtoms);
    disposeIntArray(bonded14IndexArray, num14);
    disposeRealArray(bonded14ParamArray, num14);
116
117
    if (neighborList != NULL)
        delete neighborList;
118
119
}

120
void ReferenceCalcStandardMMForceFieldKernel::initialize(const vector<vector<int> >& bondIndices, const vector<vector<double> >& bondParameters,
121
122
123
        const vector<vector<int> >& angleIndices, const vector<vector<double> >& angleParameters,
        const vector<vector<int> >& periodicTorsionIndices, const vector<vector<double> >& periodicTorsionParameters,
        const vector<vector<int> >& rbTorsionIndices, const vector<vector<double> >& rbTorsionParameters,
124
        const vector<vector<int> >& bonded14Indices, double lj14Scale, double coulomb14Scale,
125
126
        const vector<set<int> >& exclusions, const vector<vector<double> >& nonbondedParameters,
        NonbondedMethod nonbondedMethod, double nonbondedCutoff, double periodicBoxSize[3]) {
127
    numAtoms = nonbondedParameters.size();
128
129
130
131
    numBonds = bondIndices.size();
    numAngles = angleIndices.size();
    numPeriodicTorsions = periodicTorsionIndices.size();
    numRBTorsions = rbTorsionIndices.size();
132
    num14 = bonded14Indices.size();
133
134
135
136
137
138
139
140
    bondIndexArray = copyToArray(bondIndices);
    bondParamArray = copyToArray(bondParameters);
    angleIndexArray = copyToArray(angleIndices);
    angleParamArray = copyToArray(angleParameters);
    periodicTorsionIndexArray = copyToArray(periodicTorsionIndices);
    periodicTorsionParamArray = copyToArray(periodicTorsionParameters);
    rbTorsionIndexArray = copyToArray(rbTorsionIndices);
    rbTorsionParamArray = copyToArray(rbTorsionParameters);
141
    atomParamArray = allocateRealArray(numAtoms, 3);
142
    RealOpenMM sqrtEps = static_cast<RealOpenMM>( std::sqrt(138.935485) );
143
    for (int i = 0; i < numAtoms; ++i) {
144
145
146
        atomParamArray[i][0] = static_cast<RealOpenMM>( 0.5*nonbondedParameters[i][1] );
        atomParamArray[i][1] = static_cast<RealOpenMM>( 2.0*sqrt(nonbondedParameters[i][2]) );
        atomParamArray[i][2] = static_cast<RealOpenMM>( nonbondedParameters[i][0]*sqrtEps );
147
    }
148
    this->exclusions = exclusions;
149
150
151
152
153
154
155
156
157
158
159
160
161
    exclusionArray = new int*[numAtoms];
    for (int i = 0; i < numAtoms; ++i) {
        exclusionArray[i] = new int[exclusions[i].size()+1];
        exclusionArray[i][0] = exclusions[i].size();
        int index = 0;
        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
            exclusionArray[i][++index] = *iter;
    }
    bonded14IndexArray = copyToArray(bonded14Indices);
    bonded14ParamArray = allocateRealArray(num14, 3);
    for (int i = 0; i < num14; ++i) {
        int atom1 = bonded14Indices[i][0];
        int atom2 = bonded14Indices[i][1];
162
163
164
        bonded14ParamArray[i][0] = static_cast<RealOpenMM>( atomParamArray[atom1][0]+atomParamArray[atom2][0] );
        bonded14ParamArray[i][1] = static_cast<RealOpenMM>( lj14Scale*(atomParamArray[atom1][1]*atomParamArray[atom2][1]) );
        bonded14ParamArray[i][2] = static_cast<RealOpenMM>( coulomb14Scale*(atomParamArray[atom1][2]*atomParamArray[atom2][2]) );
165
    }
166
    this->nonbondedMethod = nonbondedMethod;
167
168
169
170
    this->nonbondedCutoff = (RealOpenMM) nonbondedCutoff;
    this->periodicBoxSize[0] = (RealOpenMM) periodicBoxSize[0];
    this->periodicBoxSize[1] = (RealOpenMM) periodicBoxSize[1];
    this->periodicBoxSize[2] = (RealOpenMM) periodicBoxSize[2];
171
172
173
174
175
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();
        
176
177
}

178
void ReferenceCalcStandardMMForceFieldKernel::executeForces(const Stream& positions, Stream& forces) {
179
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) positions.getImpl()).getData()); // Reference code needs to be made const correct
180
181
182
183
184
185
186
187
188
189
190
    RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) forces.getImpl()).getData();
    ReferenceBondForce refBondForce;
    ReferenceHarmonicBondIxn harmonicBond;
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, 0, 0, harmonicBond);
    ReferenceAngleBondIxn angleBond;
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, angleBond);
    ReferenceProperDihedralBond periodicTorsionBond;
    refBondForce.calculateForce(numPeriodicTorsions, periodicTorsionIndexArray, posData, periodicTorsionParamArray, forceData, 0, 0, 0, periodicTorsionBond);
    ReferenceRbDihedralBond rbTorsionBond;
    refBondForce.calculateForce(numRBTorsions, rbTorsionIndexArray, posData, rbTorsionParamArray, forceData, 0, 0, 0, rbTorsionBond);
    ReferenceLJCoulombIxn clj;
191
192
193
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
        computeNeighborListVoxelHash(*neighborList, numAtoms, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
194
        clj.setUseCutoff(nonbondedCutoff, *neighborList, 78.3f);
195
196
197
    }
    if (periodic)
        clj.setPeriodic(periodicBoxSize);
198
199
    clj.calculatePairIxn(numAtoms, posData, atomParamArray, exclusionArray, 0, forceData, 0, 0);
    ReferenceLJCoulomb14 nonbonded14;
200
    if (nonbondedMethod != NoCutoff)
201
        nonbonded14.setUseCutoff(nonbondedCutoff, 78.3f);
202
203
204
205
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
}

double ReferenceCalcStandardMMForceFieldKernel::executeEnergy(const Stream& positions) {
206
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) positions.getImpl()).getData()); // Reference code needs to be made const correct
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
    RealOpenMM** forceData = allocateRealArray(numAtoms, 3);
    int arraySize = max(max(max(max(numAtoms, numBonds), numAngles), numPeriodicTorsions), numRBTorsions);
    RealOpenMM* energyArray = new RealOpenMM[arraySize];
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceHarmonicBondIxn harmonicBond;
    for (int i = 0; i < arraySize; ++i)
        energyArray[i] = 0;
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, energyArray, 0, &energy, harmonicBond);
    ReferenceAngleBondIxn angleBond;
    for (int i = 0; i < arraySize; ++i)
        energyArray[i] = 0;
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, energyArray, 0, &energy, angleBond);
    ReferenceProperDihedralBond periodicTorsionBond;
    for (int i = 0; i < arraySize; ++i)
        energyArray[i] = 0;
    refBondForce.calculateForce(numPeriodicTorsions, periodicTorsionIndexArray, posData, periodicTorsionParamArray, forceData, energyArray, 0, &energy, periodicTorsionBond);
    ReferenceRbDihedralBond rbTorsionBond;
    for (int i = 0; i < arraySize; ++i)
        energyArray[i] = 0;
    refBondForce.calculateForce(numRBTorsions, rbTorsionIndexArray, posData, rbTorsionParamArray, forceData, energyArray, 0, &energy, rbTorsionBond);
    ReferenceLJCoulombIxn clj;
229
230
231
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
        computeNeighborListVoxelHash(*neighborList, numAtoms, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
232
        clj.setUseCutoff(nonbondedCutoff, *neighborList, 78.3f);
233
234
235
    }
    if (periodic)
        clj.setPeriodic(periodicBoxSize);
236
237
    clj.calculatePairIxn(numAtoms, posData, atomParamArray, exclusionArray, 0, forceData, 0, &energy);
    ReferenceLJCoulomb14 nonbonded14;
238
    if (nonbondedMethod != NoCutoff)
239
        nonbonded14.setUseCutoff(nonbondedCutoff, 78.3f);
240
241
242
243
244
245
246
247
    for (int i = 0; i < arraySize; ++i)
        energyArray[i] = 0;
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, energyArray, 0, &energy, nonbonded14);
    disposeRealArray(forceData, numAtoms);
    delete[] energyArray;
    return energy;
}

248
249
ReferenceCalcGBSAOBCForceFieldKernel::~ReferenceCalcGBSAOBCForceFieldKernel() {
    if (obc) {
250
        // delete obc->getObcParameters();
251
252
253
254
255
256
257
258
259
260
        delete obc;
    }
}

void ReferenceCalcGBSAOBCForceFieldKernel::initialize(const vector<vector<double> >& atomParameters, double solventDielectric, double soluteDielectric) {
    int numAtoms = atomParameters.size();
    charges.resize(numAtoms);
    vector<RealOpenMM> atomicRadii(numAtoms);
    vector<RealOpenMM> scaleFactors(numAtoms);
    for (int i = 0; i < numAtoms; ++i) {
261
262
263
        charges[i] = static_cast<RealOpenMM>( atomParameters[i][0] );
        atomicRadii[i] = static_cast<RealOpenMM>( atomParameters[i][1] );
        scaleFactors[i] = static_cast<RealOpenMM>( atomParameters[i][2] );
264
265
    }
    ObcParameters* obcParameters  = new ObcParameters(numAtoms, ObcParameters::ObcTypeII);
266
    obcParameters->setAtomicRadii(atomicRadii, SimTKOpenMMCommon::MdUnits);
267
    obcParameters->setScaledRadiusFactors(scaleFactors);
268
269
    obcParameters->setSolventDielectric( static_cast<RealOpenMM>(solventDielectric) );
    obcParameters->setSoluteDielectric( static_cast<RealOpenMM>(soluteDielectric) );
270
271
    obc = new CpuObc(obcParameters);
    obc->setIncludeAceApproximation(true);
272
273
}

274
void ReferenceCalcGBSAOBCForceFieldKernel::executeForces(const Stream& positions, Stream& forces) {
275
276
277
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) positions.getImpl()).getData()); // Reference code needs to be made const correct
    RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) forces.getImpl()).getData();
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 0);
278
279
}

280
double ReferenceCalcGBSAOBCForceFieldKernel::executeEnergy(const Stream& positions) {
281
282
283
284
285
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) positions.getImpl()).getData()); // Reference code needs to be made const correct
    RealOpenMM** forceData = allocateRealArray(positions.getSize(), 3);
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
    disposeRealArray(forceData, positions.getSize());
    return obc->getEnergy();
286
287
}

288
289
290
291
292
293
294
295
296
297
298
299
300
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
    if (dynamics)
        delete dynamics;
    if (shake)
        delete shake;
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
    if (shakeParameters)
        disposeRealArray(shakeParameters, numConstraints);
}

301
302
void ReferenceIntegrateVerletStepKernel::initialize(const vector<double>& masses, const vector<vector<int> >& constraintIndices,
        const vector<double>& constraintLengths) {
303
    this->masses = new RealOpenMM[masses.size()];
304
305
    for (size_t i = 0; i < masses.size(); ++i)
        this->masses[i] = static_cast<RealOpenMM>( masses[i] );
306
307
308
309
310
311
312
    numConstraints = constraintIndices.size();
    this->constraintIndices = allocateIntArray(numConstraints, 2);
    for (int i = 0; i < numConstraints; ++i) {
        this->constraintIndices[i][0] = constraintIndices[i][0];
        this->constraintIndices[i][1] = constraintIndices[i][1];
    }
    shakeParameters = allocateRealArray(constraintLengths.size(), 1);
313
314
    for (size_t i = 0; i < constraintLengths.size(); ++i)
        shakeParameters[i][0] = static_cast<RealOpenMM>( constraintLengths[i] );
315
316
317
}

void ReferenceIntegrateVerletStepKernel::execute(Stream& positions, Stream& velocities, const Stream& forces, double stepSize) {
318
319
320
321
322
323
324
325
326
327
    RealOpenMM** posData = ((ReferenceFloatStreamImpl&) positions.getImpl()).getData();
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) velocities.getImpl()).getData();
    RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) forces.getImpl()).getData()); // Reference code needs to be made const correct
    if (dynamics == 0 || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
            delete shake;
        }
328
        dynamics = new ReferenceVerletDynamics(positions.getSize(), static_cast<RealOpenMM>(stepSize) );
329
330
331
332
333
        shake = new ReferenceShakeAlgorithm(numConstraints, constraintIndices, shakeParameters);
        dynamics->setReferenceShakeAlgorithm(shake);
        prevStepSize = stepSize;
    }
    dynamics->update(positions.getSize(), posData, velData, forceData, masses);
334
}
335

336
337
338
339
340
341
342
343
344
345
346
347
ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
    if (shake)
        delete shake;
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
    if (shakeParameters)
        disposeRealArray(shakeParameters, numConstraints);
}
348
349
350

void ReferenceIntegrateLangevinStepKernel::initialize(const vector<double>& masses, const vector<vector<int> >& constraintIndices,
        const vector<double>& constraintLengths) {
351
    this->masses = new RealOpenMM[masses.size()];
352
353
    for (size_t i = 0; i < masses.size(); ++i)
        this->masses[i] = static_cast<RealOpenMM>( masses[i] );
354
355
356
357
358
359
360
    numConstraints = constraintIndices.size();
    this->constraintIndices = allocateIntArray(numConstraints, 2);
    for (int i = 0; i < numConstraints; ++i) {
        this->constraintIndices[i][0] = constraintIndices[i][0];
        this->constraintIndices[i][1] = constraintIndices[i][1];
    }
    shakeParameters = allocateRealArray(constraintLengths.size(), 1);
361
362
    for (size_t i = 0; i < constraintLengths.size(); ++i)
        shakeParameters[i][0] = static_cast<RealOpenMM>( constraintLengths[i] );
363
364
365
}

void ReferenceIntegrateLangevinStepKernel::execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize) {
366
367
368
369
370
371
372
373
374
375
    RealOpenMM** posData = ((ReferenceFloatStreamImpl&) positions.getImpl()).getData();
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) velocities.getImpl()).getData();
    RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) forces.getImpl()).getData()); // Reference code needs to be made const correct
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
            delete shake;
        }
376
377
378
379
380
381
        RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
        dynamics = new ReferenceStochasticDynamics(
				positions.getSize(), 
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(tau), 
				static_cast<RealOpenMM>(temperature) );
382
383
384
385
386
387
388
        shake = new ReferenceShakeAlgorithm(numConstraints, constraintIndices, shakeParameters);
        dynamics->setReferenceShakeAlgorithm(shake);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
    dynamics->update(positions.getSize(), posData, velData, forceData, masses);
389
390
}

391
392
393
394
395
396
397
398
399
400
401
402
403
ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
    if (dynamics)
        delete dynamics;
    if (shake)
        delete shake;
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
    if (shakeParameters)
        disposeRealArray(shakeParameters, numConstraints);
}

404
405
void ReferenceIntegrateBrownianStepKernel::initialize(const vector<double>& masses, const vector<vector<int> >& constraintIndices,
        const vector<double>& constraintLengths) {
406
    this->masses = new RealOpenMM[masses.size()];
407
408
    for (size_t i = 0; i < masses.size(); ++i)
        this->masses[i] = static_cast<RealOpenMM>(masses[i]);
409
410
411
412
413
414
415
    numConstraints = constraintIndices.size();
    this->constraintIndices = allocateIntArray(numConstraints, 2);
    for (int i = 0; i < numConstraints; ++i) {
        this->constraintIndices[i][0] = constraintIndices[i][0];
        this->constraintIndices[i][1] = constraintIndices[i][1];
    }
    shakeParameters = allocateRealArray(constraintLengths.size(), 1);
416
417
    for (size_t i = 0; i < constraintLengths.size(); ++i)
        shakeParameters[i][0] = static_cast<RealOpenMM>( constraintLengths[i] );
418
419
420
}

void ReferenceIntegrateBrownianStepKernel::execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize) {
421
422
423
424
425
426
427
428
429
430
    RealOpenMM** posData = ((ReferenceFloatStreamImpl&) positions.getImpl()).getData();
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) velocities.getImpl()).getData();
    RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) forces.getImpl()).getData()); // Reference code needs to be made const correct
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
            delete shake;
        }
431
432
433
434
435
        dynamics = new ReferenceBrownianDynamics(
				positions.getSize(), 
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(friction), 
				static_cast<RealOpenMM>(temperature) );
436
437
438
439
440
441
442
        shake = new ReferenceShakeAlgorithm(numConstraints, constraintIndices, shakeParameters);
        dynamics->setReferenceShakeAlgorithm(shake);
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
    dynamics->update(positions.getSize(), posData, velData, forceData, masses);
443
444
}

445
446
447
448
449
450
451
ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
    if (thermostat)
        delete thermostat;
    if (masses)
        delete[] masses;
}

452
void ReferenceApplyAndersenThermostatKernel::initialize(const vector<double>& masses) {
453
    this->masses = new RealOpenMM[masses.size()];
454
455
    for (size_t i = 0; i < masses.size(); ++i)
        this->masses[i] = static_cast<RealOpenMM>( masses[i] );
456
    thermostat = new ReferenceAndersenThermostat();
457
458
459
}

void ReferenceApplyAndersenThermostatKernel::execute(Stream& velocities, double temperature, double collisionFrequency, double stepSize) {
460
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) velocities.getImpl()).getData();
461
462
463
464
465
466
467
    thermostat->applyThermostat(
			velocities.getSize(), 
			velData, 
			masses, 
			static_cast<RealOpenMM>(temperature), 
			static_cast<RealOpenMM>(collisionFrequency), 
			static_cast<RealOpenMM>(stepSize) );
468
469
470
}

void ReferenceCalcKineticEnergyKernel::initialize(const vector<double>& masses) {
471
    this->masses = masses;
472
473
}

474
475
476
double ReferenceCalcKineticEnergyKernel::execute(const Stream& velocities) {
    RealOpenMM** velData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) velocities.getImpl()).getData()); // Reference code needs to be made const correct
    double energy = 0.0;
477
    for (size_t i = 0; i < masses.size(); ++i)
478
479
        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;
480
}
481
482
483

void ReferenceRemoveCMMotionKernel::initialize(const vector<double>& masses) {
    this->masses.resize(masses.size());
484
    for (size_t i = 0; i < masses.size(); ++i)
485
486
487
488
489
490
491
492
493
        this->masses[i] = masses[i];
}

void ReferenceRemoveCMMotionKernel::execute(Stream& velocities) {
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) velocities.getImpl()).getData();
    
    // Calculate the center of mass momentum.
    
    RealOpenMM momentum[] = {0.0, 0.0, 0.0};
494
495
496
497
    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] );
498
499
500
501
    }
    
    // Adjust the atom velocities.
    
502
503
504
505
506
507
508
    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] );
509
510
    }
}