"examples/cpp-examples/HelloArgonInFortran.f90" did not exist on "de366ef7239f73cb2dbc26357b306dfbe5edd987"
ReferenceKernels.cpp 76.3 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-2012 Stanford University and the Authors.      *
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 * Authors: Peter Eastman                                                     *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

#include "ReferenceKernels.h"
33
#include "gbsa/CpuObc.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
34
#include "gbsa/CpuGBVI.h"
35
#include "SimTKReference/ReferenceAndersenThermostat.h"
36
37
#include "SimTKReference/ReferenceAngleBondIxn.h"
#include "SimTKReference/ReferenceBondForce.h"
38
#include "SimTKReference/ReferenceBrownianDynamics.h"
39
#include "SimTKReference/ReferenceCCMAAlgorithm.h"
40
#include "SimTKReference/ReferenceCMAPTorsionIxn.h"
41
#include "SimTKReference/ReferenceCustomAngleIxn.h"
42
#include "SimTKReference/ReferenceCustomBondIxn.h"
43
#include "SimTKReference/ReferenceCustomDynamics.h"
44
#include "SimTKReference/ReferenceCustomExternalIxn.h"
45
#include "SimTKReference/ReferenceCustomGBIxn.h"
46
#include "SimTKReference/ReferenceCustomHbondIxn.h"
47
#include "SimTKReference/ReferenceCustomNonbondedIxn.h"
48
#include "SimTKReference/ReferenceCustomTorsionIxn.h"
49
50
51
#include "SimTKReference/ReferenceHarmonicBondIxn.h"
#include "SimTKReference/ReferenceLJCoulomb14.h"
#include "SimTKReference/ReferenceLJCoulombIxn.h"
52
#include "SimTKReference/ReferenceMonteCarloBarostat.h"
53
54
#include "SimTKReference/ReferenceProperDihedralBond.h"
#include "SimTKReference/ReferenceRbDihedralBond.h"
55
#include "SimTKReference/ReferenceStochasticDynamics.h"
56
57
#include "SimTKReference/ReferenceVariableStochasticDynamics.h"
#include "SimTKReference/ReferenceVariableVerletDynamics.h"
58
#include "SimTKReference/ReferenceVerletDynamics.h"
59
#include "SimTKReference/ReferenceVirtualSites.h"
60
#include "openmm/CMMotionRemover.h"
61
#include "openmm/Context.h"
62
#include "openmm/System.h"
63
#include "openmm/internal/AndersenThermostatImpl.h"
64
#include "openmm/internal/ContextImpl.h"
65
#include "openmm/internal/CustomHbondForceImpl.h"
66
#include "openmm/internal/CMAPTorsionForceImpl.h"
67
#include "openmm/internal/NonbondedForceImpl.h"
68
#include "openmm/internal/SplineFitter.h"
69
#include "openmm/Integrator.h"
70
#include "openmm/OpenMMException.h"
71
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
72
#include "lepton/CustomFunction.h"
73
74
#include "lepton/Parser.h"
#include "lepton/ParsedExpression.h"
75
#include <cmath>
76
#include <limits>
77
78
79
80

using namespace OpenMM;
using namespace std;

81
static int** allocateIntArray(int length, int width) {
82
83
84
85
86
87
    int** array = new int*[length];
    for (int i = 0; i < length; ++i)
        array[i] = new int[width];
    return array;
}

88
static RealOpenMM** allocateRealArray(int length, int width) {
89
90
91
92
93
94
    RealOpenMM** array = new RealOpenMM*[length];
    for (int i = 0; i < length; ++i)
        array[i] = new RealOpenMM[width];
    return array;
}

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

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

111
static vector<RealVec>& extractPositions(ContextImpl& context) {
112
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
113
    return *((vector<RealVec>*) data->positions);
114
115
}

116
static vector<RealVec>& extractVelocities(ContextImpl& context) {
117
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
118
    return *((vector<RealVec>*) data->velocities);
119
120
}

121
static vector<RealVec>& extractForces(ContextImpl& context) {
122
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
123
    return *((vector<RealVec>*) data->forces);
124
125
}

126
static RealVec& extractBoxSize(ContextImpl& context) {
127
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
128
    return *(RealVec*) data->periodicBoxSize;
129
130
}

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

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

148
void ReferenceCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
149
    vector<RealVec>& forceData = extractForces(context);
150
151
152
153
154
155
156
    if (includeForces) {
        int numParticles = context.getSystem().getNumParticles();
        for (int i = 0; i < numParticles; ++i) {
            forceData[i][0] = (RealOpenMM) 0.0;
            forceData[i][1] = (RealOpenMM) 0.0;
            forceData[i][2] = (RealOpenMM) 0.0;
        }
157
    }
158
159
    else
        savedForces = forceData;
160
161
}

162
double ReferenceCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
163
164
    if (!includeForces)
        extractForces(context) = savedForces; // Restore the forces so computing the energy doesn't overwrite the forces with incorrect values.
165
166
    else
        ReferenceVirtualSites::distributeForces(context.getSystem(), extractPositions(context), extractForces(context));
167
168
169
    return 0.0;
}

170
void ReferenceUpdateStateDataKernel::initialize(const System& system) {
171
172
}

173
double ReferenceUpdateStateDataKernel::getTime(const ContextImpl& context) const {
174
175
176
    return data.time;
}

177
void ReferenceUpdateStateDataKernel::setTime(ContextImpl& context, double time) {
178
179
180
    data.time = time;
}

181
182
void ReferenceUpdateStateDataKernel::getPositions(ContextImpl& context, std::vector<Vec3>& positions) {
    int numParticles = context.getSystem().getNumParticles();
183
    vector<RealVec>& posData = extractPositions(context);
184
185
186
187
188
189
190
    positions.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        positions[i] = Vec3(posData[i][0], posData[i][1], posData[i][2]);
}

void ReferenceUpdateStateDataKernel::setPositions(ContextImpl& context, const std::vector<Vec3>& positions) {
    int numParticles = context.getSystem().getNumParticles();
191
    vector<RealVec>& posData = extractPositions(context);
192
    for (int i = 0; i < numParticles; ++i) {
193
194
195
        posData[i][0] = (RealOpenMM) positions[i][0];
        posData[i][1] = (RealOpenMM) positions[i][1];
        posData[i][2] = (RealOpenMM) positions[i][2];
196
197
198
199
200
    }
}

void ReferenceUpdateStateDataKernel::getVelocities(ContextImpl& context, std::vector<Vec3>& velocities) {
    int numParticles = context.getSystem().getNumParticles();
201
    vector<RealVec>& velData = extractVelocities(context);
202
203
204
205
206
207
208
    velocities.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        velocities[i] = Vec3(velData[i][0], velData[i][1], velData[i][2]);
}

void ReferenceUpdateStateDataKernel::setVelocities(ContextImpl& context, const std::vector<Vec3>& velocities) {
    int numParticles = context.getSystem().getNumParticles();
209
    vector<RealVec>& velData = extractVelocities(context);
210
    for (int i = 0; i < numParticles; ++i) {
211
212
213
        velData[i][0] = (RealOpenMM) velocities[i][0];
        velData[i][1] = (RealOpenMM) velocities[i][1];
        velData[i][2] = (RealOpenMM) velocities[i][2];
214
215
216
217
218
    }
}

void ReferenceUpdateStateDataKernel::getForces(ContextImpl& context, std::vector<Vec3>& forces) {
    int numParticles = context.getSystem().getNumParticles();
219
    vector<RealVec>& forceData = extractForces(context);
220
221
222
223
224
    forces.resize(numParticles);
    for (int i = 0; i < numParticles; ++i)
        forces[i] = Vec3(forceData[i][0], forceData[i][1], forceData[i][2]);
}

225
void ReferenceUpdateStateDataKernel::getPeriodicBoxVectors(ContextImpl& context, Vec3& a, Vec3& b, Vec3& c) const {
226
    RealVec& box = extractBoxSize(context);
227
228
229
230
231
232
    a = Vec3(box[0], 0, 0);
    b = Vec3(0, box[1], 0);
    c = Vec3(0, 0, box[2]);
}

void ReferenceUpdateStateDataKernel::setPeriodicBoxVectors(ContextImpl& context, const Vec3& a, const Vec3& b, const Vec3& c) const {
233
    RealVec& box = extractBoxSize(context);
234
235
236
237
238
    box[0] = (RealOpenMM) a[0];
    box[1] = (RealOpenMM) b[1];
    box[2] = (RealOpenMM) c[2];
}

239
240
void ReferenceApplyConstraintsKernel::initialize(const System& system) {
    int numParticles = system.getNumParticles();
241
242
    masses.resize(numParticles);
    inverseMasses.resize(numParticles);
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
    for (int i = 0; i < numParticles; ++i) {
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
        inverseMasses[i] = 1.0/masses[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);
    }
}

ReferenceApplyConstraintsKernel::~ReferenceApplyConstraintsKernel() {
    if (constraints)
        delete constraints;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
    if (constraintDistances)
        delete[] constraintDistances;
}

void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
    if (constraints == NULL) {
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, tol);
    }
275
    vector<RealVec>& positions = extractPositions(context);
276
277
    constraints->setTolerance(tol);
    constraints->apply(data.numParticles, positions, positions, inverseMasses);
278
    ReferenceVirtualSites::computePositions(context.getSystem(), positions);
279
280
}

281
282
283
284
285
286
287
288
void ReferenceVirtualSitesKernel::initialize(const System& system) {
}

void ReferenceVirtualSitesKernel::computePositions(ContextImpl& context) {
    vector<RealVec>& positions = extractPositions(context);
    ReferenceVirtualSites::computePositions(context.getSystem(), positions);
}

289
ReferenceCalcHarmonicBondForceKernel::~ReferenceCalcHarmonicBondForceKernel() {
290
291
292
293
    disposeIntArray(bondIndexArray, numBonds);
    disposeRealArray(bondParamArray, numBonds);
}

294
void ReferenceCalcHarmonicBondForceKernel::initialize(const System& system, const HarmonicBondForce& force) {
295
296
297
298
    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
299
        int particle1, particle2;
300
        double length, k;
Peter Eastman's avatar
Peter Eastman committed
301
302
303
        force.getBondParameters(i, particle1, particle2, length, k);
        bondIndexArray[i][0] = particle1;
        bondIndexArray[i][1] = particle2;
304
305
        bondParamArray[i][0] = (RealOpenMM) length;
        bondParamArray[i][1] = (RealOpenMM) k;
306
    }
307
308
}

309
double ReferenceCalcHarmonicBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
310
311
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
312
313
314
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceHarmonicBondIxn harmonicBond;
315
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, includeEnergy ? &energy : NULL, harmonicBond);
316
317
318
    return energy;
}

319
320
321
322
323
324
325
326
327
328
329
ReferenceCalcCustomBondForceKernel::~ReferenceCalcCustomBondForceKernel() {
    disposeIntArray(bondIndexArray, numBonds);
    disposeRealArray(bondParamArray, numBonds);
}

void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const CustomBondForce& force) {
    numBonds = force.getNumBonds();
    int numParameters = force.getNumPerBondParameters();

    // Build the arrays.

330
    bondIndexArray = allocateIntArray(numBonds, 2);
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
    bondParamArray = allocateRealArray(numBonds, numParameters);
    vector<double> params;
    for (int i = 0; i < force.getNumBonds(); ++i) {
        int particle1, particle2;
        force.getBondParameters(i, particle1, particle2, params);
        bondIndexArray[i][0] = particle1;
        bondIndexArray[i][1] = particle2;
        for (int j = 0; j < numParameters; j++)
            bondParamArray[i][j] = (RealOpenMM) params[j];
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
    energyExpression = expression.createProgram();
    forceExpression = expression.differentiate("r").optimize().createProgram();
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerBondParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
}

353
double ReferenceCalcCustomBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
354
355
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
356
357
358
359
360
361
    RealOpenMM energy = 0;
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceBondForce refBondForce;
    ReferenceCustomBondIxn harmonicBond(energyExpression, forceExpression, parameterNames, globalParameters);
362
    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, includeEnergy ? &energy : NULL, harmonicBond);
363
364
365
    return energy;
}

366
367
368
369
370
371
372
373
374
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);
375
    for (int i = 0; i < force.getNumAngles(); ++i) {
Peter Eastman's avatar
Peter Eastman committed
376
        int particle1, particle2, particle3;
377
        double angle, k;
Peter Eastman's avatar
Peter Eastman committed
378
379
380
381
        force.getAngleParameters(i, particle1, particle2, particle3, angle, k);
        angleIndexArray[i][0] = particle1;
        angleIndexArray[i][1] = particle2;
        angleIndexArray[i][2] = particle3;
382
383
        angleParamArray[i][0] = (RealOpenMM) angle;
        angleParamArray[i][1] = (RealOpenMM) k;
384
    }
385
386
}

387
double ReferenceCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
388
389
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
390
391
392
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceAngleBondIxn angleBond;
393
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, includeEnergy ? &energy : NULL, angleBond);
394
395
396
    return energy;
}

397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
ReferenceCalcCustomAngleForceKernel::~ReferenceCalcCustomAngleForceKernel() {
    disposeIntArray(angleIndexArray, numAngles);
    disposeRealArray(angleParamArray, numAngles);
}

void ReferenceCalcCustomAngleForceKernel::initialize(const System& system, const CustomAngleForce& force) {
    numAngles = force.getNumAngles();
    int numParameters = force.getNumPerAngleParameters();

    // Build the arrays.

    angleIndexArray = allocateIntArray(numAngles, 3);
    angleParamArray = allocateRealArray(numAngles, numParameters);
    vector<double> params;
    for (int i = 0; i < force.getNumAngles(); ++i) {
        int particle1, particle2, particle3;
        force.getAngleParameters(i, particle1, particle2, particle3, params);
        angleIndexArray[i][0] = particle1;
        angleIndexArray[i][1] = particle2;
        angleIndexArray[i][2] = particle3;
        for (int j = 0; j < numParameters; j++)
            angleParamArray[i][j] = (RealOpenMM) params[j];
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
    energyExpression = expression.createProgram();
    forceExpression = expression.differentiate("theta").optimize().createProgram();
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerAngleParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
}

432
double ReferenceCalcCustomAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
433
434
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
435
436
437
438
439
    RealOpenMM energy = 0;
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceBondForce refBondForce;
440
    ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters);
441
    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, includeEnergy ? &energy : NULL, customAngle);
442
443
444
    return energy;
}

445
446
447
448
449
450
451
452
453
454
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
455
        int particle1, particle2, particle3, particle4, periodicity;
456
        double phase, k;
Peter Eastman's avatar
Peter Eastman committed
457
458
459
460
461
        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;
462
463
464
        torsionParamArray[i][0] = (RealOpenMM) k;
        torsionParamArray[i][1] = (RealOpenMM) phase;
        torsionParamArray[i][2] = (RealOpenMM) periodicity;
465
    }
466
467
}

468
double ReferenceCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
469
470
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
471
472
473
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceProperDihedralBond periodicTorsionBond;
474
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, periodicTorsionBond);
475
476
477
478
479
480
481
482
483
484
485
486
487
    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
488
        int particle1, particle2, particle3, particle4;
489
        double c0, c1, c2, c3, c4, c5;
Peter Eastman's avatar
Peter Eastman committed
490
491
492
493
494
        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;
495
496
497
498
499
500
        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;
501
    }
502
503
}

504
double ReferenceCalcRBTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
505
506
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
507
508
509
    RealOpenMM energy = 0;
    ReferenceBondForce refBondForce;
    ReferenceRbDihedralBond rbTorsionBond;
510
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, rbTorsionBond);
511
512
513
    return energy;
}

514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
void ReferenceCalcCMAPTorsionForceKernel::initialize(const System& system, const CMAPTorsionForce& force) {
    int numMaps = force.getNumMaps();
    int numTorsions = force.getNumTorsions();
    coeff.resize(numMaps);
    vector<double> energy;
    vector<vector<double> > c;
    for (int i = 0; i < numMaps; i++) {
        int size;
        force.getMapParameters(i, size, energy);
        CMAPTorsionForceImpl::calcMapDerivatives(size, energy, c);
        coeff[i].resize(size*size);
        for (int j = 0; j < size*size; j++) {
            coeff[i][j].resize(16);
            for (int k = 0; k < 16; k++)
                coeff[i][j][k] = c[j][k];
        }
    }
    torsionMaps.resize(numTorsions);
    torsionIndices.resize(numTorsions);
    for (int i = 0; i < numTorsions; i++) {
        torsionIndices[i].resize(8);
        force.getTorsionParameters(i, torsionMaps[i], torsionIndices[i][0], torsionIndices[i][1], torsionIndices[i][2],
            torsionIndices[i][3], torsionIndices[i][4], torsionIndices[i][5], torsionIndices[i][6], torsionIndices[i][7]);
    }
}

540
double ReferenceCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
541
542
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
543
544
545
546
547
548
    RealOpenMM totalEnergy = 0;
    ReferenceCMAPTorsionIxn torsion(coeff, torsionMaps, torsionIndices);
    torsion.calculateIxn(posData, forceData, &totalEnergy);
    return totalEnergy;
}

549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
ReferenceCalcCustomTorsionForceKernel::~ReferenceCalcCustomTorsionForceKernel() {
    disposeIntArray(torsionIndexArray, numTorsions);
    disposeRealArray(torsionParamArray, numTorsions);
}

void ReferenceCalcCustomTorsionForceKernel::initialize(const System& system, const CustomTorsionForce& force) {
    numTorsions = force.getNumTorsions();
    int numParameters = force.getNumPerTorsionParameters();

    // Build the arrays.

    torsionIndexArray = allocateIntArray(numTorsions, 4);
    torsionParamArray = allocateRealArray(numTorsions, numParameters);
    vector<double> params;
    for (int i = 0; i < force.getNumTorsions(); ++i) {
        int particle1, particle2, particle3, particle4;
        force.getTorsionParameters(i, particle1, particle2, particle3, particle4, params);
        torsionIndexArray[i][0] = particle1;
        torsionIndexArray[i][1] = particle2;
        torsionIndexArray[i][2] = particle3;
        torsionIndexArray[i][3] = particle4;
        for (int j = 0; j < numParameters; j++)
            torsionParamArray[i][j] = (RealOpenMM) params[j];
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
    energyExpression = expression.createProgram();
    forceExpression = expression.differentiate("theta").optimize().createProgram();
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerTorsionParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
}

585
double ReferenceCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
586
587
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
588
589
590
591
592
593
    RealOpenMM energy = 0;
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceBondForce refBondForce;
    ReferenceCustomTorsionIxn customTorsion(energyExpression, forceExpression, parameterNames, globalParameters);
594
    refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, includeEnergy ? &energy : NULL, customTorsion);
595
596
597
    return energy;
}

598
ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
Peter Eastman's avatar
Peter Eastman committed
599
600
    disposeRealArray(particleParamArray, numParticles);
    disposeIntArray(exclusionArray, numParticles);
601
602
603
604
605
606
    disposeIntArray(bonded14IndexArray, num14);
    disposeRealArray(bonded14ParamArray, num14);
    if (neighborList != NULL)
        delete neighborList;
}

607
608
609
610
void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {

    // Identify which exceptions are 1-4 interactions.

Peter Eastman's avatar
Peter Eastman committed
611
    numParticles = force.getNumParticles();
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
    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();
627
628
    bonded14IndexArray = allocateIntArray(num14, 2);
    bonded14ParamArray = allocateRealArray(num14, 3);
Peter Eastman's avatar
Peter Eastman committed
629
630
    particleParamArray = allocateRealArray(numParticles, 3);
    for (int i = 0; i < numParticles; ++i) {
631
        double charge, radius, depth;
Peter Eastman's avatar
Peter Eastman committed
632
633
634
        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));
635
        particleParamArray[i][2] = static_cast<RealOpenMM>(charge);
636
    }
637
    this->exclusions = exclusions;
Peter Eastman's avatar
Peter Eastman committed
638
639
    exclusionArray = new int*[numParticles];
    for (int i = 0; i < numParticles; ++i) {
640
641
642
643
644
645
646
        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
647
        int particle1, particle2;
648
        double charge, radius, depth;
649
        force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth);
Peter Eastman's avatar
Peter Eastman committed
650
651
        bonded14IndexArray[i][0] = particle1;
        bonded14IndexArray[i][1] = particle2;
652
653
        bonded14ParamArray[i][0] = static_cast<RealOpenMM>(radius);
        bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
654
        bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge);
655
    }
656
    nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
657
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
658
659
660
661
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();
662
663
664
    if (nonbondedMethod == Ewald) {
        double alpha;
        NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmax[0], kmax[1], kmax[2]);
665
        ewaldAlpha = (RealOpenMM) alpha;
666
667
668
669
    }
    else if (nonbondedMethod == PME) {
        double alpha;
        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2]);
670
        ewaldAlpha = (RealOpenMM) alpha;
671
    }
672
    rfDielectric = (RealOpenMM)force.getReactionFieldDielectric();
673
674
675
676
    if (force.getUseDispersionCorrection())
        dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
    else
        dispersionCoefficient = 0.0;
677
678
}

679
double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
680
681
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
682
683
    RealOpenMM energy = 0;
    ReferenceLJCoulombIxn clj;
684
    bool periodic = (nonbondedMethod == CutoffPeriodic);
685
    bool ewald  = (nonbondedMethod == Ewald);
686
    bool pme  = (nonbondedMethod == PME);
687
    if (nonbondedMethod != NoCutoff) {
688
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic || ewald || pme, nonbondedCutoff, 0.0);
689
        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
690
    }
691
    if (periodic || ewald || pme)
692
        clj.setPeriodic(extractBoxSize(context));
693
694
    if (ewald)
        clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
695
    if (pme)
696
        clj.setUsePME(ewaldAlpha, gridSize);
697
698
699
700
701
702
    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, includeEnergy ? &energy : NULL, includeDirect, includeReciprocal);
    if (includeDirect) {
        ReferenceBondForce refBondForce;
        ReferenceLJCoulomb14 nonbonded14;
        refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
    }
703
    if (periodic || ewald || pme) {
704
        RealVec& boxSize = extractBoxSize(context);
705
706
        energy += dispersionCoefficient/(boxSize[0]*boxSize[1]*boxSize[2]);
    }
707
708
709
    return energy;
}

710
class ReferenceTabulatedFunction : public Lepton::CustomFunction {
711
public:
712
713
714
715
716
717
718
    ReferenceTabulatedFunction(double min, double max, const vector<double>& values) :
            min(min), max(max), values(values) {
        int numValues = values.size();
        x.resize(numValues);
        for (int i = 0; i < numValues; i++)
            x[i] = min+i*(max-min)/(numValues-1);
        SplineFitter::createNaturalSpline(x, values, derivs);
719
720
721
722
723
    }
    int getNumArguments() const {
        return 1;
    }
    double evaluate(const double* arguments) const {
724
725
        double t = arguments[0];
        if (t < min || t > max)
726
            return 0.0;
727
        return SplineFitter::evaluateSpline(x, values, derivs, t);
728
729
    }
    double evaluateDerivative(const double* arguments, const int* derivOrder) const {
730
731
        double t = arguments[0];
        if (t < min || t > max)
732
            return 0.0;
733
        return SplineFitter::evaluateSplineDerivative(x, values, derivs, t);
734
735
    }
    CustomFunction* clone() const {
736
        return new ReferenceTabulatedFunction(min, max, values);
737
738
    }
    double min, max;
739
    vector<double> x, values, derivs;
740
741
};

742
743
744
745
746
747
748
749
750
ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
    disposeRealArray(particleParamArray, numParticles);
    disposeIntArray(exclusionArray, numParticles);
    if (neighborList != NULL)
        delete neighborList;
}

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

751
    // Record the exclusions.
752
753
754

    numParticles = force.getNumParticles();
    exclusions.resize(numParticles);
755
    for (int i = 0; i < force.getNumExclusions(); i++) {
756
        int particle1, particle2;
757
        force.getExclusionParticles(i, particle1, particle2);
758
759
760
761
762
763
        exclusions[particle1].insert(particle2);
        exclusions[particle2].insert(particle1);
    }

    // Build the arrays.

764
    int numParameters = force.getNumPerParticleParameters();
765
766
    particleParamArray = allocateRealArray(numParticles, numParameters);
    for (int i = 0; i < numParticles; ++i) {
767
        vector<double> parameters;
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
        force.getParticleParameters(i, parameters);
        for (int j = 0; j < numParameters; j++)
            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
    }
    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;
    }
    nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();

787
788
789
790
791
792
793
    // 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;
794
795
        force.getFunctionParameters(i, name, values, min, max);
        functions[name] = new ReferenceTabulatedFunction(min, max, values);
796
797
    }

798
799
    // Parse the various expressions used to calculate the force.

800
    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
801
802
    energyExpression = expression.createProgram();
    forceExpression = expression.differentiate("r").optimize().createProgram();
803
804
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerParticleParameterName(i));
805
806
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
807
808
809
810
811

    // Delete the custom functions.

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

814
double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
815
816
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
817
    RealOpenMM energy = 0;
818
    ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames);
819
820
    bool periodic = (nonbondedMethod == CutoffPeriodic);
    if (nonbondedMethod != NoCutoff) {
821
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic, nonbondedCutoff, 0.0);
822
823
824
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    if (periodic)
825
        ixn.setPeriodic(extractBoxSize(context));
826
    map<string, double> globalParameters;
827
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
828
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
829
    ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, includeEnergy ? &energy : NULL);
830
831
832
    return energy;
}

833
ReferenceCalcGBSAOBCForceKernel::~ReferenceCalcGBSAOBCForceKernel() {
834
    if (obc) {
835
        // delete obc->getObcParameters();
836
837
838
839
        delete obc;
    }
}

840
void ReferenceCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCForce& force) {
Peter Eastman's avatar
Peter Eastman committed
841
842
843
844
845
    int numParticles = system.getNumParticles();
    charges.resize(numParticles);
    vector<RealOpenMM> atomicRadii(numParticles);
    vector<RealOpenMM> scaleFactors(numParticles);
    for (int i = 0; i < numParticles; ++i) {
846
        double charge, radius, scalingFactor;
Peter Eastman's avatar
Peter Eastman committed
847
        force.getParticleParameters(i, charge, radius, scalingFactor);
848
849
850
        charges[i] = static_cast<RealOpenMM>(charge);
        atomicRadii[i] = static_cast<RealOpenMM>(radius);
        scaleFactors[i] = static_cast<RealOpenMM>(scalingFactor);
851
    }
852
    ObcParameters* obcParameters = new ObcParameters(numParticles, ObcParameters::ObcTypeII);
853
    obcParameters->setAtomicRadii(atomicRadii);
854
    obcParameters->setScaledRadiusFactors(scaleFactors);
855
856
    obcParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
    obcParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
857
858
    if (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff)
        obcParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
859
    isPeriodic = (force.getNonbondedMethod() == GBSAOBCForce::CutoffPeriodic);
860
861
    obc = new CpuObc(obcParameters);
    obc->setIncludeAceApproximation(true);
862
863
}

864
double ReferenceCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
865
866
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
867
868
    if (isPeriodic)
        obc->getObcParameters()->setPeriodic(extractBoxSize(context));
Mark Friedrichs's avatar
Mark Friedrichs committed
869
    return obc->computeBornEnergyForces(posData, charges, forceData);
870
871
}

Mark Friedrichs's avatar
Mark Friedrichs committed
872
873
ReferenceCalcGBVIForceKernel::~ReferenceCalcGBVIForceKernel() {
    if (gbvi) {
Mark Friedrichs's avatar
Mark Friedrichs committed
874
875
        GBVIParameters * gBVIParameters = gbvi->getGBVIParameters();
        delete gBVIParameters;
Mark Friedrichs's avatar
Mark Friedrichs committed
876
877
878
879
880
        delete gbvi;
    }
}

void ReferenceCalcGBVIForceKernel::initialize(const System& system, const GBVIForce& force, const std::vector<double> & inputScaledRadii ) {
881

Mark Friedrichs's avatar
Mark Friedrichs committed
882
    int numParticles = system.getNumParticles();
883

Mark Friedrichs's avatar
Mark Friedrichs committed
884
885
886
887
    charges.resize(numParticles);
    vector<RealOpenMM> atomicRadii(numParticles);
    vector<RealOpenMM> scaledRadii(numParticles);
    vector<RealOpenMM> gammas(numParticles);
888

Mark Friedrichs's avatar
Mark Friedrichs committed
889
890
891
892
893
894
895
896
    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]);
    }
897

Mark Friedrichs's avatar
Mark Friedrichs committed
898
    GBVIParameters * gBVIParameters = new GBVIParameters(numParticles);
899

Mark Friedrichs's avatar
Mark Friedrichs committed
900
901
902
    gBVIParameters->setAtomicRadii(atomicRadii);
    gBVIParameters->setGammaParameters(gammas);
    gBVIParameters->setScaledRadii(scaledRadii);
903
904
    gBVIParameters->setSolventDielectric(static_cast<RealOpenMM>(force.getSolventDielectric()));
    gBVIParameters->setSoluteDielectric(static_cast<RealOpenMM>(force.getSoluteDielectric()));
905

906
907
908
    gBVIParameters->setBornRadiusScalingMethod(force.getBornRadiusScalingMethod());
    gBVIParameters->setQuinticUpperBornRadiusLimit(static_cast<RealOpenMM>(force.getQuinticUpperBornRadiusLimit()));
    gBVIParameters->setQuinticLowerLimitFactor(static_cast<RealOpenMM>(force.getQuinticLowerLimitFactor()));
909

910
911
    if (force.getNonbondedMethod() != GBVIForce::NoCutoff)
        gBVIParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
912
    isPeriodic = (force.getNonbondedMethod() == GBVIForce::CutoffPeriodic);
Mark Friedrichs's avatar
Mark Friedrichs committed
913
914
915
    gbvi = new CpuGBVI(gBVIParameters);
}

916
double ReferenceCalcGBVIForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
Mark Friedrichs's avatar
Mark Friedrichs committed
917

918
    vector<RealVec>& posData = extractPositions(context);
Mark Friedrichs's avatar
Mark Friedrichs committed
919

920
921
    if (isPeriodic)
        gbvi->getGBVIParameters()->setPeriodic(extractBoxSize(context));
Mark Friedrichs's avatar
Mark Friedrichs committed
922
923

    RealOpenMM energy;
924
    if (includeForces) {
925
        vector<RealVec>& forceData = extractForces(context);
Mark Friedrichs's avatar
Mark Friedrichs committed
926
927
928
929
930
        gbvi->computeBornForces(posData, charges, forceData);
        energy = 0.0;
    }
    if( includeEnergy ){
        energy = gbvi->computeBornEnergy(posData, charges);
931
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
932
933
934
    return static_cast<double>(energy);
}

935
936
937
938
939
940
941
ReferenceCalcCustomGBForceKernel::~ReferenceCalcCustomGBForceKernel() {
    disposeRealArray(particleParamArray, numParticles);
    if (neighborList != NULL)
        delete neighborList;
}

void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
942
943
944
945
946
947
948
949
950
951
952
953
    if (force.getNumComputedValues() > 0) {
        string name, expression;
        CustomGBForce::ComputationType type;
        force.getComputedValueParameters(0, name, expression, type);
        if (type == CustomGBForce::SingleParticle)
            throw OpenMMException("ReferencePlatform requires that the first computed value for a CustomGBForce be of type ParticlePair or ParticlePairNoExclusions.");
        for (int i = 1; i < force.getNumComputedValues(); i++) {
            force.getComputedValueParameters(i, name, expression, type);
            if (type != CustomGBForce::SingleParticle)
                throw OpenMMException("ReferencePlatform requires that a CustomGBForce only have one computed value of type ParticlePair or ParticlePairNoExclusions.");
        }
    }
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989

    // Record the exclusions.

    numParticles = force.getNumParticles();
    exclusions.resize(numParticles);
    for (int i = 0; i < force.getNumExclusions(); i++) {
        int particle1, particle2;
        force.getExclusionParticles(i, particle1, particle2);
        exclusions[particle1].insert(particle2);
        exclusions[particle2].insert(particle1);
    }

    // Build the arrays.

    int numPerParticleParameters = force.getNumPerParticleParameters();
    particleParamArray = allocateRealArray(numParticles, numPerParticleParameters);
    for (int i = 0; i < numParticles; ++i) {
        vector<double> parameters;
        force.getParticleParameters(i, parameters);
        for (int j = 0; j < numPerParticleParameters; j++)
            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
    }
    for (int i = 0; i < numPerParticleParameters; i++)
        particleParameterNames.push_back(force.getPerParticleParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
    nonbondedMethod = CalcCustomGBForceKernel::NonbondedMethod(force.getNonbondedMethod());
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
    if (nonbondedMethod == NoCutoff)
        neighborList = NULL;
    else
        neighborList = new NeighborList();

    // Create custom functions for the tabulated functions.

    map<string, Lepton::CustomFunction*> functions;
990
991
992
993
    for (int i = 0; i < force.getNumFunctions(); i++) {
        string name;
        vector<double> values;
        double min, max;
994
995
        force.getFunctionParameters(i, name, values, min, max);
        functions[name] = new ReferenceTabulatedFunction(min, max, values);
996
    }
997
998
999

    // Parse the expressions for computed values.

1000
    valueDerivExpressions.resize(force.getNumComputedValues());
1001
    valueGradientExpressions.resize(force.getNumComputedValues());
1002
1003
1004
1005
1006
1007
1008
1009
    for (int i = 0; i < force.getNumComputedValues(); i++) {
        string name, expression;
        CustomGBForce::ComputationType type;
        force.getComputedValueParameters(i, name, expression, type);
        Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
        valueExpressions.push_back(ex.createProgram());
        valueTypes.push_back(type);
        valueNames.push_back(name);
1010
1011
1012
        if (i == 0)
            valueDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
        else {
1013
1014
1015
            valueGradientExpressions[i].push_back(ex.differentiate("x").optimize().createProgram());
            valueGradientExpressions[i].push_back(ex.differentiate("y").optimize().createProgram());
            valueGradientExpressions[i].push_back(ex.differentiate("z").optimize().createProgram());
1016
1017
1018
            for (int j = 0; j < i; j++)
                valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
        }
1019
1020
    }

1021
    // Parse the expressions for energy terms.
1022
1023

    energyDerivExpressions.resize(force.getNumEnergyTerms());
1024
    energyGradientExpressions.resize(force.getNumEnergyTerms());
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
    for (int i = 0; i < force.getNumEnergyTerms(); i++) {
        string expression;
        CustomGBForce::ComputationType type;
        force.getEnergyTermParameters(i, expression, type);
        Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
        energyExpressions.push_back(ex.createProgram());
        energyTypes.push_back(type);
        if (type != CustomGBForce::SingleParticle)
            energyDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
        for (int j = 0; j < force.getNumComputedValues(); j++) {
1035
            if (type == CustomGBForce::SingleParticle) {
1036
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
1037
1038
1039
1040
                energyGradientExpressions[i].push_back(ex.differentiate("x").optimize().createProgram());
                energyGradientExpressions[i].push_back(ex.differentiate("y").optimize().createProgram());
                energyGradientExpressions[i].push_back(ex.differentiate("z").optimize().createProgram());
            }
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
            else {
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").optimize().createProgram());
                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").optimize().createProgram());
            }
        }
    }

    // Delete the custom functions.

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

1054
double ReferenceCalcCustomGBForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
1055
1056
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
1057
    RealOpenMM energy = 0;
1058
1059
    ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
        energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
1060
    bool periodic = (nonbondedMethod == CutoffPeriodic);
1061
1062
    if (periodic)
        ixn.setPeriodic(extractBoxSize(context));
1063
    if (nonbondedMethod != NoCutoff) {
1064
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic, nonbondedCutoff, 0.0);
1065
1066
1067
1068
1069
        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
    }
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
1070
    ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, includeEnergy ? &energy : NULL);
1071
1072
1073
    return energy;
}

1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
ReferenceCalcCustomExternalForceKernel::~ReferenceCalcCustomExternalForceKernel() {
    disposeRealArray(particleParamArray, numParticles);
}

void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, const CustomExternalForce& force) {
    numParticles = force.getNumParticles();
    int numParameters = force.getNumPerParticleParameters();

    // Build the arrays.

    particles.resize(numParticles);
    particleParamArray = allocateRealArray(numParticles, numParameters);
    vector<double> params;
    for (int i = 0; i < numParticles; ++i) {
        force.getParticleParameters(i, particles[i], params);
        for (int j = 0; j < numParameters; j++)
            particleParamArray[i][j] = (RealOpenMM) params[j];
    }

    // Parse the expression used to calculate the force.

    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
    energyExpression = expression.createProgram();
    forceExpressionX = expression.differentiate("x").optimize().createProgram();
    forceExpressionY = expression.differentiate("y").optimize().createProgram();
    forceExpressionZ = expression.differentiate("z").optimize().createProgram();
    for (int i = 0; i < numParameters; i++)
        parameterNames.push_back(force.getPerParticleParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
}

1106
double ReferenceCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
1107
1108
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
1109
1110
1111
1112
1113
1114
    RealOpenMM energy = 0;
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
    ReferenceCustomExternalIxn force(energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ, parameterNames, globalParameters);
    for (int i = 0; i < numParticles; ++i)
1115
        force.calculateForce(particles[i], posData, particleParamArray[i], forceData, includeEnergy ? &energy : NULL);
1116
1117
1118
    return energy;
}

1119
1120
1121
1122
ReferenceCalcCustomHbondForceKernel::~ReferenceCalcCustomHbondForceKernel() {
    disposeRealArray(donorParamArray, numDonors);
    disposeRealArray(acceptorParamArray, numAcceptors);
    disposeIntArray(exclusionArray, numDonors);
1123
1124
    if (ixn != NULL)
        delete ixn;
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
}

void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const CustomHbondForce& force) {

    // Record the exclusions.

    numDonors = force.getNumDonors();
    numAcceptors = force.getNumAcceptors();
    numParticles = system.getNumParticles();
    exclusions.resize(numDonors);
    for (int i = 0; i < force.getNumExclusions(); i++) {
        int donor, acceptor;
        force.getExclusionParticles(i, donor, acceptor);
        exclusions[donor].insert(acceptor);
    }

    // Build the arrays.

1143
    vector<vector<int> > donorParticles(numDonors);
1144
1145
1146
1147
    int numDonorParameters = force.getNumPerDonorParameters();
    donorParamArray = allocateRealArray(numDonors, numDonorParameters);
    for (int i = 0; i < numDonors; ++i) {
        vector<double> parameters;
1148
1149
1150
1151
1152
        int d1, d2, d3;
        force.getDonorParameters(i, d1, d2, d3, parameters);
        donorParticles[i].push_back(d1);
        donorParticles[i].push_back(d2);
        donorParticles[i].push_back(d3);
1153
1154
1155
        for (int j = 0; j < numDonorParameters; j++)
            donorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
    }
1156
    vector<vector<int> > acceptorParticles(numAcceptors);
1157
1158
1159
1160
    int numAcceptorParameters = force.getNumPerAcceptorParameters();
    acceptorParamArray = allocateRealArray(numAcceptors, numAcceptorParameters);
    for (int i = 0; i < numAcceptors; ++i) {
        vector<double> parameters;
1161
1162
1163
1164
1165
        int a1, a2, a3;
        force.getAcceptorParameters(i, a1, a2, a3, parameters);
        acceptorParticles[i].push_back(a1);
        acceptorParticles[i].push_back(a2);
        acceptorParticles[i].push_back(a3);
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
        for (int j = 0; j < numAcceptorParameters; j++)
            acceptorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
    }
    exclusionArray = new int*[numDonors];
    for (int i = 0; i < numDonors; ++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;
    }
1177
    NonbondedMethod nonbondedMethod = CalcCustomHbondForceKernel::NonbondedMethod(force.getNonbondedMethod());
1178
1179
1180
1181
1182
1183
1184
1185
1186
    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();

    // 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;
1187
1188
        force.getFunctionParameters(i, name, values, min, max);
        functions[name] = new ReferenceTabulatedFunction(min, max, values);
1189
1190
    }

1191
    // Parse the expression and create the object used to calculate the interaction.
1192

1193
1194
1195
    map<string, vector<int> > distances;
    map<string, vector<int> > angles;
    map<string, vector<int> > dihedrals;
1196
    Lepton::ParsedExpression energyExpression = CustomHbondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
1197
1198
    vector<string> donorParameterNames;
    vector<string> acceptorParameterNames;
1199
1200
1201
1202
1203
1204
    for (int i = 0; i < numDonorParameters; i++)
        donorParameterNames.push_back(force.getPerDonorParameterName(i));
    for (int i = 0; i < numAcceptorParameters; i++)
        acceptorParameterNames.push_back(force.getPerAcceptorParameterName(i));
    for (int i = 0; i < force.getNumGlobalParameters(); i++)
        globalParameterNames.push_back(force.getGlobalParameterName(i));
1205
    ixn = new ReferenceCustomHbondIxn(donorParticles, acceptorParticles, energyExpression, donorParameterNames, acceptorParameterNames, distances, angles, dihedrals);
1206
    isPeriodic = (nonbondedMethod == CutoffPeriodic);
1207
1208
    if (nonbondedMethod != NoCutoff)
        ixn->setUseCutoff(nonbondedCutoff);
1209
1210
1211
1212
1213
1214
1215

    // Delete the custom functions.

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

1216
double ReferenceCalcCustomHbondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
1217
1218
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& forceData = extractForces(context);
1219
1220
    if (isPeriodic)
        ixn->setPeriodic(extractBoxSize(context));
1221
1222
1223
1224
    RealOpenMM energy = 0;
    map<string, double> globalParameters;
    for (int i = 0; i < (int) globalParameterNames.size(); i++)
        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
1225
    ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
1226
1227
1228
    return energy;
}

1229
1230
1231
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
    if (dynamics)
        delete dynamics;
1232
1233
    if (constraints)
        delete constraints;
1234
1235
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
1236
1237
    if (constraintDistances)
        delete[] constraintDistances;
1238
1239
}

1240
void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
1241
    int numParticles = system.getNumParticles();
1242
    masses.resize(numParticles);
1243
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1244
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1245
1246
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
1247
    constraintDistances = new RealOpenMM[numConstraints];
1248
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
1249
        int particle1, particle2;
1250
        double distance;
Peter Eastman's avatar
Peter Eastman committed
1251
1252
1253
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
1254
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
1255
    }
1256
1257
}

1258
void ReferenceIntegrateVerletStepKernel::execute(ContextImpl& context, const VerletIntegrator& integrator) {
1259
    double stepSize = integrator.getStepSize();
1260
1261
1262
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& velData = extractVelocities(context);
    vector<RealVec>& forceData = extractForces(context);
1263
1264
1265
1266
1267
    if (dynamics == 0 || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
1268
            delete constraints;
1269
        }
Peter Eastman's avatar
Peter Eastman committed
1270
        dynamics = new ReferenceVerletDynamics(context.getSystem().getNumParticles(), static_cast<RealOpenMM>(stepSize) );
1271
1272
1273
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
1274
        dynamics->setReferenceConstraintAlgorithm(constraints);
1275
1276
        prevStepSize = stepSize;
    }
1277
    constraints->setTolerance(integrator.getConstraintTolerance());
1278
    dynamics->update(context.getSystem(), posData, velData, forceData, masses);
1279
    data.time += stepSize;
1280
    data.stepCount++;
1281
}
1282

1283
1284
1285
ReferenceIntegrateLangevinStepKernel::~ReferenceIntegrateLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
1286
1287
    if (constraints)
        delete constraints;
1288
1289
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
1290
1291
    if (constraintDistances)
        delete[] constraintDistances;
1292
}
1293

1294
void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
1295
    int numParticles = system.getNumParticles();
1296
    masses.resize(numParticles);
1297
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1298
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1299
1300
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
1301
    constraintDistances = new RealOpenMM[numConstraints];
1302
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
1303
        int particle1, particle2;
1304
        double distance;
Peter Eastman's avatar
Peter Eastman committed
1305
1306
1307
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
1308
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
1309
    }
1310
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
1311
1312
}

1313
void ReferenceIntegrateLangevinStepKernel::execute(ContextImpl& context, const LangevinIntegrator& integrator) {
1314
1315
1316
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
1317
1318
1319
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& velData = extractVelocities(context);
    vector<RealVec>& forceData = extractForces(context);
1320
1321
1322
1323
1324
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
1325
            delete constraints;
1326
        }
1327
1328
        RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
        dynamics = new ReferenceStochasticDynamics(
Peter Eastman's avatar
Peter Eastman committed
1329
				context.getSystem().getNumParticles(), 
1330
1331
1332
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(tau), 
				static_cast<RealOpenMM>(temperature) );
1333
1334
1335
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
1336
        dynamics->setReferenceConstraintAlgorithm(constraints);
1337
1338
1339
1340
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
1341
    constraints->setTolerance(integrator.getConstraintTolerance());
1342
    dynamics->update(context.getSystem(), posData, velData, forceData, masses);
1343
    data.time += stepSize;
1344
    data.stepCount++;
1345
1346
}

1347
1348
1349
ReferenceIntegrateBrownianStepKernel::~ReferenceIntegrateBrownianStepKernel() {
    if (dynamics)
        delete dynamics;
1350
1351
    if (constraints)
        delete constraints;
1352
1353
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
1354
1355
    if (constraintDistances)
        delete[] constraintDistances;
1356
1357
}

1358
void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
Peter Eastman's avatar
Peter Eastman committed
1359
    int numParticles = system.getNumParticles();
1360
    masses.resize(numParticles);
1361
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1362
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1363
1364
    numConstraints = system.getNumConstraints();
    constraintIndices = allocateIntArray(numConstraints, 2);
1365
    constraintDistances = new RealOpenMM[numConstraints];
1366
    for (int i = 0; i < numConstraints; ++i) {
Peter Eastman's avatar
Peter Eastman committed
1367
        int particle1, particle2;
1368
        double distance;
Peter Eastman's avatar
Peter Eastman committed
1369
1370
1371
        system.getConstraintParameters(i, particle1, particle2, distance);
        constraintIndices[i][0] = particle1;
        constraintIndices[i][1] = particle2;
1372
        constraintDistances[i] = static_cast<RealOpenMM>(distance);
1373
    }
1374
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
1375
1376
}

1377
void ReferenceIntegrateBrownianStepKernel::execute(ContextImpl& context, const BrownianIntegrator& integrator) {
1378
1379
1380
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double stepSize = integrator.getStepSize();
1381
1382
1383
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& velData = extractVelocities(context);
    vector<RealVec>& forceData = extractForces(context);
1384
1385
1386
1387
1388
    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
        // Recreate the computation objects with the new parameters.
        
        if (dynamics) {
            delete dynamics;
1389
            delete constraints;
1390
        }
1391
        dynamics = new ReferenceBrownianDynamics(
Peter Eastman's avatar
Peter Eastman committed
1392
				context.getSystem().getNumParticles(), 
1393
1394
1395
				static_cast<RealOpenMM>(stepSize), 
				static_cast<RealOpenMM>(friction), 
				static_cast<RealOpenMM>(temperature) );
1396
1397
1398
        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
        findAnglesForCCMA(context.getSystem(), angles);
        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
1399
        dynamics->setReferenceConstraintAlgorithm(constraints);
1400
1401
1402
1403
        prevTemp = temperature;
        prevFriction = friction;
        prevStepSize = stepSize;
    }
1404
    constraints->setTolerance(integrator.getConstraintTolerance());
1405
    dynamics->update(context.getSystem(), posData, velData, forceData, masses);
1406
    data.time += stepSize;
1407
    data.stepCount++;
1408
1409
}

1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
ReferenceIntegrateVariableLangevinStepKernel::~ReferenceIntegrateVariableLangevinStepKernel() {
    if (dynamics)
        delete dynamics;
    if (constraints)
        delete constraints;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
    if (constraintDistances)
        delete[] constraintDistances;
}

void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& system, const VariableLangevinIntegrator& integrator) {
    int numParticles = system.getNumParticles();
1423
    masses.resize(numParticles);
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
    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());
}

1440
void ReferenceIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
1441
1442
1443
    double temperature = integrator.getTemperature();
    double friction = integrator.getFriction();
    double errorTol = integrator.getErrorTolerance();
1444
1445
1446
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& velData = extractVelocities(context);
    vector<RealVec>& forceData = extractForces(context);
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
    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;
    }
1464
    constraints->setTolerance(integrator.getConstraintTolerance());
1465
    RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
1466
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize);
1467
1468
1469
1470
1471
1472
    data.time += dynamics->getDeltaT();
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
}

1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
ReferenceIntegrateVariableVerletStepKernel::~ReferenceIntegrateVariableVerletStepKernel() {
    if (dynamics)
        delete dynamics;
    if (constraints)
        delete constraints;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
    if (constraintDistances)
        delete[] constraintDistances;
}

void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system, const VariableVerletIntegrator& integrator) {
    int numParticles = system.getNumParticles();
1486
    masses.resize(numParticles);
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
    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);
    }
}

1502
void ReferenceIntegrateVariableVerletStepKernel::execute(ContextImpl& context, const VariableVerletIntegrator& integrator, double maxTime) {
1503
    double errorTol = integrator.getErrorTolerance();
1504
1505
1506
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& velData = extractVelocities(context);
    vector<RealVec>& forceData = extractForces(context);
1507
    if (dynamics == 0 || errorTol != prevErrorTol) {
1508
1509
1510
1511
1512
1513
        // Recreate the computation objects with the new parameters.

        if (dynamics) {
            delete dynamics;
            delete constraints;
        }
1514
        dynamics = new ReferenceVariableVerletDynamics(context.getSystem().getNumParticles(), (RealOpenMM) errorTol);
1515
1516
1517
1518
        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);
1519
        prevErrorTol = errorTol;
1520
    }
1521
    constraints->setTolerance(integrator.getConstraintTolerance());
1522
    RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
1523
    dynamics->update(context.getSystem(), posData, velData, forceData, masses, maxStepSize);
1524
    data.time += dynamics->getDeltaT();
1525
1526
1527
    if (dynamics->getDeltaT() == maxStepSize)
        data.time = maxTime; // Avoid round-off error
    data.stepCount++;
1528
1529
}

1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
ReferenceIntegrateCustomStepKernel::~ReferenceIntegrateCustomStepKernel() {
    if (dynamics)
        delete dynamics;
    if (constraints)
        delete constraints;
    if (constraintIndices)
        disposeIntArray(constraintIndices, numConstraints);
    if (constraintDistances)
        delete[] constraintDistances;
}

void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const CustomIntegrator& integrator) {
    int numParticles = system.getNumParticles();
    masses.resize(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);
    }
    perDofValues.resize(integrator.getNumPerDofVariables());
    for (int i = 0; i < (int) perDofValues.size(); i++)
        perDofValues[i].resize(numParticles);
}

void ReferenceIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegrator& integrator, bool& forcesAreValid) {
    vector<RealVec>& posData = extractPositions(context);
    vector<RealVec>& velData = extractVelocities(context);
    vector<RealVec>& forceData = extractForces(context);
    if (dynamics == 0) {
        // Create the computation objects.

        dynamics = new ReferenceCustomDynamics(context.getSystem().getNumParticles(), integrator);
        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);
    }
    
    // Record global variables.
    
    map<string, double> globals;
    globals["dt"] = integrator.getStepSize();
    for (int i = 0; i < integrator.getNumGlobalVariables(); i++)
        globals[integrator.getGlobalVariableName(i)] = globalValues[i];
    
    // Execute the step.
    
    constraints->setTolerance(integrator.getConstraintTolerance());
    dynamics->update(context, context.getSystem().getNumParticles(), posData, velData, forceData, masses, globals, perDofValues, forcesAreValid);
    
    // Record changed global variables.
    
    integrator.setStepSize(globals["dt"]);
    for (int i = 0; i < (int) globalValues.size(); i++)
        globalValues[i] = globals[integrator.getGlobalVariableName(i)];
    data.time += dynamics->getDeltaT();
    data.stepCount++;
}

void ReferenceIntegrateCustomStepKernel::getGlobalVariables(ContextImpl& context, vector<double>& values) const {
    values = globalValues;
}

void ReferenceIntegrateCustomStepKernel::setGlobalVariables(ContextImpl& context, const vector<double>& values) {
    globalValues = values;
}

void ReferenceIntegrateCustomStepKernel::getPerDofVariable(ContextImpl& context, int variable, vector<Vec3>& values) const {
    values.resize(perDofValues[variable].size());
    for (int i = 0; i < (int) values.size(); i++)
        values[i] = perDofValues[variable][i];
}

void ReferenceIntegrateCustomStepKernel::setPerDofVariable(ContextImpl& context, int variable, const vector<Vec3>& values) {
    perDofValues[variable].resize(values.size());
    for (int i = 0; i < (int) values.size(); i++)
        perDofValues[variable][i] = values[i];
}

1617
1618
1619
1620
1621
ReferenceApplyAndersenThermostatKernel::~ReferenceApplyAndersenThermostatKernel() {
    if (thermostat)
        delete thermostat;
}

1622
void ReferenceApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
Peter Eastman's avatar
Peter Eastman committed
1623
    int numParticles = system.getNumParticles();
1624
    masses.resize(numParticles);
1625
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1626
        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
1627
    this->thermostat = new ReferenceAndersenThermostat();
1628
    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
1629
    particleGroups = AndersenThermostatImpl::calcParticleGroups(system);
1630
1631
}

1632
void ReferenceApplyAndersenThermostatKernel::execute(ContextImpl& context) {
1633
    vector<RealVec>& velData = extractVelocities(context);
1634
1635
1636
1637
    thermostat->applyThermostat(particleGroups, velData, masses,
        static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::Temperature())),
        static_cast<RealOpenMM>(context.getParameter(AndersenThermostat::CollisionFrequency())),
        static_cast<RealOpenMM>(context.getIntegrator().getStepSize()));
1638
1639
}

1640
1641
1642
1643
1644
1645
1646
1647
1648
ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel() {
    if (barostat)
        delete barostat;
}

void ReferenceApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& barostat) {
}

void ReferenceApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
1649
1650
    if (barostat == NULL)
        barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), context.getMolecules());
1651
    vector<RealVec>& posData = extractPositions(context);
1652
    RealVec& boxSize = extractBoxSize(context);
1653
1654
1655
1656
    barostat->applyBarostat(posData, boxSize, scale);
}

void ReferenceApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
1657
    vector<RealVec>& posData = extractPositions(context);
1658
1659
1660
    barostat->restorePositions(posData);
}

1661
void ReferenceCalcKineticEnergyKernel::initialize(const System& system) {
Peter Eastman's avatar
Peter Eastman committed
1662
1663
    int numParticles = system.getNumParticles();
    masses.resize(numParticles);
1664
    for (int i = 0; i < numParticles; ++i)
Peter Eastman's avatar
Peter Eastman committed
1665
        masses[i] = system.getParticleMass(i);
1666
1667
}

1668
double ReferenceCalcKineticEnergyKernel::execute(ContextImpl& context) {
1669
    vector<RealVec>& velData = extractVelocities(context);
1670
    double energy = 0.0;
1671
    for (size_t i = 0; i < masses.size(); ++i)
1672
1673
        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;
1674
}
1675

1676
1677
void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
    frequency = force.getFrequency();
Peter Eastman's avatar
Peter Eastman committed
1678
    masses.resize(system.getNumParticles());
1679
    for (size_t i = 0; i < masses.size(); ++i)
Peter Eastman's avatar
Peter Eastman committed
1680
        masses[i] = system.getParticleMass(i);
1681
1682
}

1683
void ReferenceRemoveCMMotionKernel::execute(ContextImpl& context) {
1684
    if (data.stepCount%frequency != 0)
1685
        return;
1686
    vector<RealVec>& velData = extractVelocities(context);
1687
1688
1689
1690
    
    // Calculate the center of mass momentum.
    
    RealOpenMM momentum[] = {0.0, 0.0, 0.0};
1691
    RealOpenMM mass = 0.0;
1692
    for (size_t i = 0; i < masses.size(); ++i) {
1693
1694
1695
1696
        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]);
        mass += static_cast<RealOpenMM>(masses[i]);
1697
1698
    }
    
Peter Eastman's avatar
Peter Eastman committed
1699
    // Adjust the particle velocities.
1700
    
1701
1702
1703
    momentum[0] /= mass;
    momentum[1] /= mass;
    momentum[2] /= mass;
1704
    for (size_t i = 0; i < masses.size(); ++i) {
1705
1706
1707
1708
1709
        if (masses[i] != 0.0) {
            velData[i][0] -= momentum[0];
            velData[i][1] -= momentum[1];
            velData[i][2] -= momentum[2];
        }
1710
1711
    }
}