ReferenceKernels.cpp 55.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/ReferenceAndersenThermostat.h"
37
38
#include "SimTKReference/ReferenceAngleBondIxn.h"
#include "SimTKReference/ReferenceBondForce.h"
39
#include "SimTKReference/ReferenceBrownianDynamics.h"
40
#include "SimTKReference/ReferenceCCMAAlgorithm.h"
41
#include "SimTKReference/ReferenceCustomNonbondedIxn.h"
42
43
44
45
46
#include "SimTKReference/ReferenceHarmonicBondIxn.h"
#include "SimTKReference/ReferenceLJCoulomb14.h"
#include "SimTKReference/ReferenceLJCoulombIxn.h"
#include "SimTKReference/ReferenceProperDihedralBond.h"
#include "SimTKReference/ReferenceRbDihedralBond.h"
47
#include "SimTKReference/ReferenceStochasticDynamics.h"
48
49
#include "SimTKReference/ReferenceVariableStochasticDynamics.h"
#include "SimTKReference/ReferenceVariableVerletDynamics.h"
50
#include "SimTKReference/ReferenceVerletDynamics.h"
51
52
#include "openmm/CMMotionRemover.h"
#include "openmm/System.h"
53
#include "openmm/internal/ContextImpl.h"
54
#include "openmm/Integrator.h"
55
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
56
#include "lepton/CustomFunction.h"
57
58
#include "lepton/Parser.h"
#include "lepton/ParsedExpression.h"
59
#include <cmath>
60
#include <limits>
61
62
63
64

using namespace OpenMM;
using namespace std;

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

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

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

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

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

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

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

129
130
131
void ReferenceInitializeForcesKernel::initialize(const System& system) {
}

132
void ReferenceInitializeForcesKernel::execute(ContextImpl& context) {
133
134
135
136
    double zero[] = {0.0, 0.0, 0.0};
    context.getForces().fillWithValue(zero);
}

137
138
139
void ReferenceUpdateTimeKernel::initialize(const System& system) {
}

140
double ReferenceUpdateTimeKernel::getTime(const ContextImpl& context) const {
141
142
143
    return data.time;
}

144
void ReferenceUpdateTimeKernel::setTime(ContextImpl& context, double time) {
145
146
147
    data.time = time;
}

148
ReferenceCalcHarmonicBondForceKernel::~ReferenceCalcHarmonicBondForceKernel() {
149
150
151
152
    disposeIntArray(bondIndexArray, numBonds);
    disposeRealArray(bondParamArray, numBonds);
}

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

168
void ReferenceCalcHarmonicBondForceKernel::executeForces(ContextImpl& context) {
169
170
171
172
173
174
175
    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);
}

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

212
void ReferenceCalcHarmonicAngleForceKernel::executeForces(ContextImpl& context) {
213
214
215
216
217
218
219
    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);
}

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

258
void ReferenceCalcPeriodicTorsionForceKernel::executeForces(ContextImpl& context) {
259
260
261
262
263
264
265
    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);
}

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

307
void ReferenceCalcRBTorsionForceKernel::executeForces(ContextImpl& context) {
308
309
310
311
312
313
314
    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);
}

315
double ReferenceCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
316
    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
317
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
318
319
320
321
322
323
324
    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
325
    disposeRealArray(forceData, context.getSystem().getNumParticles());
326
327
328
329
330
    delete[] energyArray;
    return energy;
}

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

339
340
341
342
void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {

    // Identify which exceptions are 1-4 interactions.

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

427
void ReferenceCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
428
429
    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();
430
    ReferenceLJCoulombIxn clj;
431
    bool periodic = (nonbondedMethod == CutoffPeriodic);
432
    bool ewald  = (nonbondedMethod == Ewald);
433
    bool pme  = (nonbondedMethod == PME);
434
    if (nonbondedMethod != NoCutoff) {
435
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
436
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
437
    }
438
    if (periodic||ewald||pme)
439
        clj.setPeriodic(periodicBoxSize);
440
441
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
442
    if (pme)
443
        clj.setUsePME(ewaldAlpha, gridSize);
444
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
445
    ReferenceBondForce refBondForce;
446
    ReferenceLJCoulomb14 nonbonded14;
Peter Eastman's avatar
Peter Eastman committed
447
    if (nonbondedMethod == CutoffNonPeriodic || nonbondedMethod == CutoffPeriodic)
448
        nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
449
450
451
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
}

452
double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
453
    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
454
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
455
456
    RealOpenMM energy = 0;
    ReferenceLJCoulombIxn clj;
457
    bool periodic = (nonbondedMethod == CutoffPeriodic);
458
    bool ewald  = (nonbondedMethod == Ewald);
459
    bool pme  = (nonbondedMethod == PME);
460
    if (nonbondedMethod != NoCutoff) {
461
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
462
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
463
    }
464
    if (periodic || ewald || pme)
465
        clj.setPeriodic(periodicBoxSize);
466
467
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
468
    if (pme)
469
        clj.setUsePME(ewaldAlpha, gridSize);
470
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
471
    ReferenceBondForce refBondForce;
472
    ReferenceLJCoulomb14 nonbonded14;
Peter Eastman's avatar
Peter Eastman committed
473
    if (nonbondedMethod == CutoffNonPeriodic || nonbondedMethod == CutoffPeriodic)
474
        nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
Peter Eastman's avatar
Peter Eastman committed
475
476
    RealOpenMM* energyArray = new RealOpenMM[num14];
    for (int i = 0; i < num14; ++i)
477
478
        energyArray[i] = 0;
    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, energyArray, 0, &energy, nonbonded14);
Peter Eastman's avatar
Peter Eastman committed
479
    disposeRealArray(forceData, numParticles);
480
481
482
483
    delete[] energyArray;
    return energy;
}

484
485
486
487
488
489
490
491
492
493
494
495
496
class ReferenceCalcCustomNonbondedForceKernel::TabulatedFunction : public Lepton::CustomFunction {
public:
    TabulatedFunction(double min, double max, const vector<double>& values, bool interpolating) :
            min(min), max(max), values(values), interpolating(interpolating) {
    }
    int getNumArguments() const {
        return 1;
    }
    /**
     * Given the function argument, find the local spline coefficients.
     */
    void findCoefficients(double& x, double* coeff) const {
        int length = values.size();
497
498
        double scale = (length-1)/(max-min);
        int index = std::floor((x-min)*scale);
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
        double points[4];
        points[0] = (index == 0 ? 2*values[0]-values[1] : values[index-1]);
        points[1] = values[index];
        points[2] = (index > length-2 ? values[length-1] : values[index+1]);
        points[3] = (index > length-3 ? 2*values[length-1]-values[length-2] : values[index+2]);
        if (interpolating) {
            coeff[0] = points[1];
            coeff[1] = 0.5*(-points[0]+points[2]);
            coeff[2] = 0.5*(2.0*points[0]-5.0*points[1]+4.0*points[2]-points[3]);
            coeff[3] = 0.5*(-points[0]+3.0*points[1]-3.0*points[2]+points[3]);
        }
        else {
            coeff[0] = (points[0]+4.0*points[1]+points[2])/6.0;
            coeff[1] = (-3.0*points[0]+3.0*points[2])/6.0;
            coeff[2] = (3.0*points[0]-6.0*points[1]+3.0*points[2])/6.0;
            coeff[3] = (-points[0]+3.0*points[1]-3.0*points[2]+points[3])/6.0;
        }
516
        x = (x-min)*scale-index;
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
    }
    double evaluate(const double* arguments) const {
        double x = arguments[0];
        if (x < min || x > max)
            return 0.0;
        double coeff[4];
        findCoefficients(x, coeff);
        return coeff[0]+x*(coeff[1]+x*(coeff[2]+x*coeff[3]));
    }
    double evaluateDerivative(const double* arguments, const int* derivOrder) const {
        double x = arguments[0];
        if (x < min || x > max)
            return 0.0;
        double coeff[4];
        findCoefficients(x, coeff);
        double scale = (values.size()-1)/(max-min);
533
        return scale*(coeff[1]+x*(2.0*coeff[2]+x*3.0*coeff[3])); // We assume a first derivative, because that's the only order ever used by CustomNonbondedForce.
534
535
536
537
538
539
540
541
542
    }
    CustomFunction* clone() const {
        return new TabulatedFunction(min, max, values, interpolating);
    }
    double min, max;
    vector<double> values;
    bool interpolating;
};

543
544
545
ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
    disposeRealArray(particleParamArray, numParticles);
    disposeIntArray(exclusionArray, numParticles);
546
547
    disposeIntArray(exceptionIndexArray, numExceptions);
    disposeRealArray(exceptionParamArray, numExceptions);
548
549
550
551
552
553
    if (neighborList != NULL)
        delete neighborList;
}

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

554
    // Identify which exceptions are actual interactions.
555
556
557

    numParticles = force.getNumParticles();
    exclusions.resize(numParticles);
558
    vector<int> exceptions;
559
560
561
562
563
564
565
    vector<double> parameters;
    for (int i = 0; i < force.getNumExceptions(); i++) {
        int particle1, particle2;
        force.getExceptionParameters(i, particle1, particle2, parameters);
        exclusions[particle1].insert(particle2);
        exclusions[particle2].insert(particle1);
        if (parameters.size() > 0)
566
            exceptions.push_back(i);
567
568
569
570
    }

    // Build the arrays.

571
    numExceptions = exceptions.size();
572
    int numParameters = force.getNumParameters();
573
574
    exceptionIndexArray = allocateIntArray(numExceptions, 2);
    exceptionParamArray = allocateRealArray(numExceptions, numParameters);
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
    particleParamArray = allocateRealArray(numParticles, numParameters);
    for (int i = 0; i < numParticles; ++i) {
        force.getParticleParameters(i, parameters);
        for (int j = 0; j < numParameters; j++)
            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
    }
    this->exclusions = exclusions;
    exclusionArray = new int*[numParticles];
    for (int i = 0; i < numParticles; ++i) {
        exclusionArray[i] = new int[exclusions[i].size()+1];
        exclusionArray[i][0] = exclusions[i].size();
        int index = 0;
        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
            exclusionArray[i][++index] = *iter;
    }
590
    for (int i = 0; i < numExceptions; ++i) {
591
        int particle1, particle2;
592
593
594
        force.getExceptionParameters(exceptions[i], particle1, particle2, parameters);
        exceptionIndexArray[i][0] = particle1;
        exceptionIndexArray[i][1] = particle2;
595
        for (int j = 0; j < numParameters; j++)
596
            exceptionParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
597
598
599
600
    }
    nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
    Vec3 boxVectors[3];
601
    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
602
603
604
605
606
607
608
609
    periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
    periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
    periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();

610
611
612
613
614
615
616
617
618
619
620
621
    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
    for (int i = 0; i < force.getNumFunctions(); i++) {
        string name;
        vector<double> values;
        double min, max;
        bool interpolating;
        force.getFunctionParameters(i, name, values, min, max, interpolating);
        functions[name] = new TabulatedFunction(min, max, values, interpolating);
    }

622
623
    // Parse the various expressions used to calculate the force.

624
    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
625
626
    energyExpression = expression.createProgram();
    forceExpression = expression.differentiate("r").optimize().createProgram();
627
628
    for (int i = 0; i < numParameters; i++) {
        parameterNames.push_back(force.getParameterName(i));
629
        combiningRules.push_back(Lepton::Parser::parse(force.getParameterCombiningRule(i), functions).optimize().createProgram());
630
631
632
    }
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
633
634
635
636
637

    // Delete the custom functions.

    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
        delete iter->second;
638
639
640
641
642
}

void ReferenceCalcCustomNonbondedForceKernel::executeForces(ContextImpl& 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();
643
    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, combiningRules);
644
645
646
647
648
649
650
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    if (periodic)
        ixn.setPeriodic(periodicBoxSize);
651
    map<string, double> globalParameters;
652
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
653
654
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, 0);
655
    ixn.calculateExceptionIxn(numExceptions, exceptionIndexArray, posData, exceptionParamArray, globalParameters, forceData, 0, 0);
656
657
658
659
660
661
}

double ReferenceCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData()); // Reference code needs to be made const correct
    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
    RealOpenMM energy = 0;
662
    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, combiningRules);
663
664
665
666
667
668
669
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    if (periodic)
        ixn.setPeriodic(periodicBoxSize);
670
    map<string, double> globalParameters;
671
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
672
673
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, &energy);
674
    ixn.calculateExceptionIxn(numExceptions, exceptionIndexArray, posData, exceptionParamArray, globalParameters, forceData, 0, &energy);
675
676
677
678
    disposeRealArray(forceData, numParticles);
    return energy;
}

679
ReferenceCalcGBSAOBCForceKernel::~ReferenceCalcGBSAOBCForceKernel() {
680
    if (obc) {
681
        // delete obc->getObcParameters();
682
683
684
685
        delete obc;
    }
}

686
void ReferenceCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
Peter Eastman's avatar
Peter Eastman committed
687
688
689
690
691
    int numParticles = system.getNumParticles();
    charges.resize(numParticles);
    vector<RealOpenMM> atomicRadii(numParticles);
    vector<RealOpenMM> scaleFactors(numParticles);
    for (int i = 0; i < numParticles; ++i) {
692
        double charge, radius, scalingFactor;
Peter Eastman's avatar
Peter Eastman committed
693
        force.getParticleParameters(i, charge, radius, scalingFactor);
694
695
696
        charges[i] = static_cast<RealOpenMM>(charge);
        atomicRadii[i] = static_cast<RealOpenMM>(radius);
        scaleFactors[i] = static_cast<RealOpenMM>(scalingFactor);
697
    }
Peter Eastman's avatar
Peter Eastman committed
698
    ObcParameters* obcParameters  = new ObcParameters(numParticles, ObcParameters::ObcTypeII);
699
    obcParameters->setAtomicRadii(atomicRadii);
700
    obcParameters->setScaledRadiusFactors(scaleFactors);
701
702
    obcParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
    obcParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
703
704
705
706
707
708
709
710
711
712
    if (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff)
        obcParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
    if (force.getNonbondedMethod() == GBSAOBCForce::CutoffPeriodic) {
        Vec3 boxVectors[3];
        system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
        RealOpenMM periodicBoxSize[3];
        periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
        periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
        periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
        obcParameters->setPeriodic(periodicBoxSize);
713
    }
714
715
    obc = new CpuObc(obcParameters);
    obc->setIncludeAceApproximation(true);
716
717
}

718
void ReferenceCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
719
720
    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();
721
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
722
723
}

724
double ReferenceCalcGBSAOBCForceKernel::executeEnergy(ContextImpl& context) {
725
    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
726
    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
727
    obc->computeImplicitSolventForces(posData, &charges[0], forceData, 1);
Peter Eastman's avatar
Peter Eastman committed
728
    disposeRealArray(forceData, context.getSystem().getNumParticles());
729
    return obc->getEnergy();
730
731
}

Mark Friedrichs's avatar
Mark Friedrichs committed
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
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);
756
757
758
759
760
761
762
763
764
765
766
767
    gBVIParameters->setSolventDielectric(static_cast<RealOpenMM>(force.getSolventDielectric()));
    gBVIParameters->setSoluteDielectric(static_cast<RealOpenMM>(force.getSoluteDielectric()));
    if (force.getNonbondedMethod() != GBVIForce::NoCutoff)
        gBVIParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
    if (force.getNonbondedMethod() == GBVIForce::CutoffPeriodic) {
        Vec3 boxVectors[3];
        system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
        RealOpenMM periodicBoxSize[3];
        periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
        periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
        periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
        gBVIParameters->setPeriodic(periodicBoxSize);
Mark Friedrichs's avatar
Mark Friedrichs committed
768
769
770
771
    }
    gbvi = new CpuGBVI(gBVIParameters);
}

772
void ReferenceCalcGBVIForceKernel::executeForces(ContextImpl& context) {
Mark Friedrichs's avatar
Mark Friedrichs committed
773
774
775
776
777
778
779
780
781

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

782
double ReferenceCalcGBVIForceKernel::executeEnergy(ContextImpl& context) {
Mark Friedrichs's avatar
Mark Friedrichs committed
783
784
785
786
787
788
789
790
    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);
}

791
792
793
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
    if (dynamics)
        delete dynamics;
794
795
    if (constraints)
        delete constraints;
796
797
798
799
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
800
801
    if (constraintDistances)
        delete[] constraintDistances;
802
803
}

804
void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
805
806
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
807
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
808
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
809
810
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
811
    constraintDistances = new RealOpenMM[numConstraints];
812
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
813
        int particle1, particle2;
814
        double distance;
Peter Eastman's avatar
Peter Eastman committed
815
816
817
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
818
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
819
    }
820
821
}

822
void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
823
824
825
826
    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
827
828
829
830
831
    if (dynamics == 0 || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
832
            delete constraints;
833
        }
Peter Eastman's avatar
Peter Eastman committed
834
        dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize) );
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());
838
        dynamics->setReferenceConstraintAlgorithm(constraints);
839
840
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
841
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
842
    data.time += stepSize;
843
    data.stepCount++;
844
}
845

846
847
848
ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
849
850
    if (constraints)
        delete constraints;
851
852
853
854
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
855
856
    if (constraintDistances)
        delete[] constraintDistances;
857
}
858

859
void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
860
861
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
862
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
863
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
864
865
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
866
    constraintDistances = new RealOpenMM[numConstraints];
867
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
868
        int particle1, particle2;
869
        double distance;
Peter Eastman's avatar
Peter Eastman committed
870
871
872
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
873
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
874
    }
875
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
876
877
}

878
void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
879
880
881
882
883
884
    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
885
886
887
888
889
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
890
            delete constraints;
891
        }
892
893
        RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
        dynamics = new ReferenceStochasticDynamics(
Peter Eastman's avatar
Peter Eastman committed
894
				context.getSystem().getNumParticles(), 
895
896
897
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(tau), 
				static_cast<RealOpenMM>(temperature) );
898
899
900
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
901
        dynamics->setReferenceConstraintAlgorithm(constraints);
902
903
904
905
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
906
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
907
    data.time += stepSize;
908
    data.stepCount++;
909
910
}

911
912
913
ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
    if (dynamics)
        delete dynamics;
914
915
    if (constraints)
        delete constraints;
916
917
918
919
    if (masses)
        delete[] masses;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
920
921
    if (constraintDistances)
        delete[] constraintDistances;
922
923
}

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

943
void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
944
945
946
947
948
949
    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
950
951
952
953
954
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
955
            delete constraints;
956
        }
957
        dynamics = new ReferenceBrownianDynamics(
Peter Eastman's avatar
Peter Eastman committed
958
				context.getSystem().getNumParticles(), 
959
960
961
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(friction), 
				static_cast<RealOpenMM>(temperature) );
962
963
964
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
965
        dynamics->setReferenceConstraintAlgorithm(constraints);
966
967
968
969
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
Peter Eastman's avatar
Peter Eastman committed
970
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
971
    data.time += stepSize;
972
    data.stepCount++;
973
974
}

975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
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());
}

1007
void ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
    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++;
}

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

1070
void ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
1071
    double errorTol = integrator.getErrorTolerance();
1072
1073
1074
    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
1075
    if (dynamics == 0 || errorTol != prevErrorTol) {
1076
1077
1078
1079
1080
1081
        // Recreate the computation objects with the new parameters.

        if (dynamics) {
            delete dynamics;
            delete constraints;
        }
1082
        dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), (RealOpenMM) errorTol);
1083
1084
1085
1086
        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);
1087
        prevErrorTol = errorTol;
1088
    }
1089
1090
    RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
1091
    data.time += dynamics->getDeltaT();
1092
1093
1094
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
1095
1096
}

1097
1098
1099
1100
1101
1102
1103
ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
    if (thermostat)
        delete thermostat;
    if (masses)
        delete[] masses;
}

1104
void ReferenceApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
Peter Eastman's avatar
Peter Eastman committed
1105
1106
    int numParticles = system.getNumParticles();
    masses = new RealOpenMM[numParticles];
1107
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1108
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1109
    this->thermostat = new ReferenceAndersenThermostat();
1110
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
1111
1112
}

1113
void ReferenceApplyAndersenThermostatKernel::execute(ContextImpl& context) {
1114
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
1115
    thermostat->applyThermostat(
1116
			context.getVelocities().getSize(), 
1117
1118
			velData, 
			masses, 
1119
1120
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::Temperature())), 
			static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::CollisionFrequency())), 
1121
			static_cast<RealOpenMM>(context.getIntegrator().getStepSize()) );
1122
1123
}

1124
void ReferenceCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
1125
1126
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
1127
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1128
        masses[i] = system.getParticleMass(i);
1129
1130
}

1131
double ReferenceCalcKineticEnergyKernel::execute(ContextImpl& context) {
1132
    RealOpenMM** velData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData()); // Reference code needs to be made const correct
1133
    double energy = 0.0;
1134
    for (size_t i = 0; i < masses.size(); ++i)
1135
1136
        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;
1137
}
1138

1139
1140
void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
    frequency = force.getFrequency();
Peter Eastman's avatar
Peter Eastman committed
1141
    masses.resize(system.getNumParticles());
1142
    for (size_t i = 0; i < masses.size(); ++i)
Peter Eastman's avatar
Peter Eastman committed
1143
        masses[i] = system.getParticleMass(i);
1144
1145
}

1146
void ReferenceRemoveCMMotionKernel::execute(ContextImpl& context) {
1147
    if (data.stepCount%frequency != 0)
1148
1149
        return;
    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
1150
1151
1152
1153
    
    // Calculate the center of mass momentum.
    
    RealOpenMM momentum[] = {0.0, 0.0, 0.0};
1154
    RealOpenMM mass = 0.0;
1155
1156
1157
1158
    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] );
1159
        mass += static_cast<RealOpenMM>( masses[i] );
1160
1161
    }
    
Peter Eastman's avatar
Peter Eastman committed
1162
    // Adjust the particle velocities.
1163
    
1164
1165
1166
    momentum[0] /= mass;
    momentum[1] /= mass;
    momentum[2] /= mass;
1167
    for (size_t i = 0; i < masses.size(); ++i) {
1168
1169
1170
        velData[i][0] -= momentum[0];
        velData[i][1] -= momentum[1];
        velData[i][2] -= momentum[2];
1171
1172
    }
}