AmoebaReferenceKernels.cpp 41.3 KB
Newer Older
1
/* -------------------------------------------------------------------------- *
2
 *                               OpenMMAmoeba                                 *
3
4
5
6
7
8
 * -------------------------------------------------------------------------- *
 * 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-2020 Stanford University and the Authors.      *
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
 * Authors:                                                                   *
 * Contributors:                                                              *
 *                                                                            *
 * This program is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU Lesser General Public License as published   *
 * by the Free Software Foundation, either version 3 of the License, or       *
 * (at your option) any later version.                                        *
 *                                                                            *
 * This program is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
 * GNU Lesser General Public License for more details.                        *
 *                                                                            *
 * You should have received a copy of the GNU Lesser General Public License   *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
 * -------------------------------------------------------------------------- */

#include "AmoebaReferenceKernels.h"
28
#include "AmoebaReferenceTorsionTorsionForce.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
29
#include "AmoebaReferenceWcaDispersionForce.h"
30
#include "AmoebaReferenceGeneralizedKirkwoodForce.h"
31
#include "openmm/internal/AmoebaTorsionTorsionForceImpl.h"
32
#include "openmm/internal/AmoebaWcaDispersionForceImpl.h"
33
34
#include "ReferencePlatform.h"
#include "openmm/internal/ContextImpl.h"
35
#include "openmm/AmoebaMultipoleForce.h"
peastman's avatar
peastman committed
36
#include "openmm/HippoNonbondedForce.h"
37
#include "openmm/internal/AmoebaMultipoleForceImpl.h"
38
#include "openmm/internal/AmoebaVdwForceImpl.h"
39
#include "openmm/internal/AmoebaGeneralizedKirkwoodForceImpl.h"
40
41
#include "openmm/NonbondedForce.h"
#include "openmm/internal/NonbondedForceImpl.h"
peastman's avatar
peastman committed
42
#include "SimTKReference/AmoebaReferenceHippoNonbondedForce.h"
43
44
45
46
47
48
49
50
51

#include <cmath>
#ifdef _MSC_VER
#include <windows.h>
#endif

using namespace OpenMM;
using namespace std;

peastman's avatar
peastman committed
52
static vector<Vec3>& extractPositions(ContextImpl& context) {
53
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
54
    return *data->positions;
55
}
56

peastman's avatar
peastman committed
57
static vector<Vec3>& extractVelocities(ContextImpl& context) {
58
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
59
    return *data->velocities;
60
}
61

peastman's avatar
peastman committed
62
static vector<Vec3>& extractForces(ContextImpl& context) {
63
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
64
    return *data->forces;
65
}
66

peastman's avatar
peastman committed
67
static Vec3& extractBoxSize(ContextImpl& context) {
68
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
69
    return *data->periodicBoxSize;
70
71
}

peastman's avatar
peastman committed
72
static Vec3* extractBoxVectors(ContextImpl& context) {
73
    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
74
    return data->periodicBoxVectors;
75
76
}

77
78
// ***************************************************************************

79
ReferenceCalcAmoebaTorsionTorsionForceKernel::ReferenceCalcAmoebaTorsionTorsionForceKernel(const std::string& name, const Platform& platform, const System& system) :
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
                CalcAmoebaTorsionTorsionForceKernel(name, platform), system(system) {
}

ReferenceCalcAmoebaTorsionTorsionForceKernel::~ReferenceCalcAmoebaTorsionTorsionForceKernel() {
}

void ReferenceCalcAmoebaTorsionTorsionForceKernel::initialize(const System& system, const AmoebaTorsionTorsionForce& force) {

    numTorsionTorsions = force.getNumTorsionTorsions();

    // torsion-torsion parameters

    for (int ii = 0; ii < numTorsionTorsions; ii++) {
        int particle1Index, particle2Index, particle3Index, particle4Index, particle5Index, chiralCheckAtomIndex, gridIndex;
        force.getTorsionTorsionParameters(ii, particle1Index, particle2Index, particle3Index,
                                          particle4Index, particle5Index, chiralCheckAtomIndex, gridIndex);
96
97
98
99
100
101
102
        particle1.push_back(particle1Index); 
        particle2.push_back(particle2Index); 
        particle3.push_back(particle3Index); 
        particle4.push_back(particle4Index); 
        particle5.push_back(particle5Index); 
        chiralCheckAtom.push_back(chiralCheckAtomIndex); 
        gridIndices.push_back(gridIndex); 
103
    }
104
    usePeriodic = force.usesPeriodicBoundaryConditions();
105
106
107
108
109
110
111

    // torsion-torsion grids

    numTorsionTorsionGrids = force.getNumTorsionTorsionGrids();
    torsionTorsionGrids.resize(numTorsionTorsionGrids);
    for (int ii = 0; ii < numTorsionTorsionGrids; ii++) {

112
113
        const TorsionTorsionGrid grid = force.getTorsionTorsionGrid(ii);
        torsionTorsionGrids[ii].resize(grid.size());
114
115
116
117
118

        // check if grid needs to be reordered: x-angle should be 'slow' index

        TorsionTorsionGrid reorderedGrid;
        int reorder = 0; 
119
120
        if (grid[0][0][0] != grid[0][1][0]) {
            AmoebaTorsionTorsionForceImpl::reorderGrid(grid, reorderedGrid);
121
122
123
            reorder = 1; 
        }    

124
125
        for (unsigned int kk = 0; kk < grid.size(); kk++) {

126
            torsionTorsionGrids[ii][kk].resize(grid[kk].size());
127
128
            for (unsigned int jj = 0; jj < grid[kk].size(); jj++) {

129
130
                torsionTorsionGrids[ii][kk][jj].resize(grid[kk][jj].size());
                if (reorder) {
131
                    for (unsigned int ll = 0; ll < grid[ll][jj].size(); ll++) {
peastman's avatar
peastman committed
132
                        torsionTorsionGrids[ii][kk][jj][ll] = reorderedGrid[kk][jj][ll];
133
134
135
                    }
                } else {
                    for (unsigned int ll = 0; ll < grid[ll][jj].size(); ll++) {
peastman's avatar
peastman committed
136
                        torsionTorsionGrids[ii][kk][jj][ll] = grid[kk][jj][ll];
137
                    }
138
139
140
141
142
143
144
145
                }
            }
        }
    }
}

double ReferenceCalcAmoebaTorsionTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {

peastman's avatar
peastman committed
146
147
    vector<Vec3>& posData   = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
148
    AmoebaReferenceTorsionTorsionForce amoebaReferenceTorsionTorsionForce;
149
150
    if (usePeriodic)
        amoebaReferenceTorsionTorsionForce.setPeriodic(extractBoxVectors(context));
peastman's avatar
peastman committed
151
    double energy = amoebaReferenceTorsionTorsionForce.calculateForceAndEnergy(numTorsionTorsions, posData,
152
                                                                                         particle1, particle2, particle3, particle4, particle5,
153
                                                                                         chiralCheckAtom, gridIndices, torsionTorsionGrids, forceData);
154
155
156
    return static_cast<double>(energy);
}

157
158
159
160
/* -------------------------------------------------------------------------- *
 *                             AmoebaMultipole                                *
 * -------------------------------------------------------------------------- */

161
ReferenceCalcAmoebaMultipoleForceKernel::ReferenceCalcAmoebaMultipoleForceKernel(const std::string& name, const Platform& platform, const System& system) :
162
163
         CalcAmoebaMultipoleForceKernel(name, platform), system(system), numMultipoles(0), mutualInducedMaxIterations(60), mutualInducedTargetEpsilon(1.0e-03),
                                                         usePme(false),alphaEwald(0.0), cutoffDistance(1.0) {  
164

165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
}

ReferenceCalcAmoebaMultipoleForceKernel::~ReferenceCalcAmoebaMultipoleForceKernel() {
}

void ReferenceCalcAmoebaMultipoleForceKernel::initialize(const System& system, const AmoebaMultipoleForce& force) {

    numMultipoles   = force.getNumMultipoles();

    charges.resize(numMultipoles);
    dipoles.resize(3*numMultipoles);
    quadrupoles.resize(9*numMultipoles);
    tholes.resize(numMultipoles);
    dampingFactors.resize(numMultipoles);
    polarity.resize(numMultipoles);
    axisTypes.resize(numMultipoles);
181
182
183
    multipoleAtomZs.resize(numMultipoles);
    multipoleAtomXs.resize(numMultipoles);
    multipoleAtomYs.resize(numMultipoles);
184
185
186
187
188
    multipoleAtomCovalentInfo.resize(numMultipoles);

    int dipoleIndex      = 0;
    int quadrupoleIndex  = 0;
    double totalCharge   = 0.0;
189
    for (int ii = 0; ii < numMultipoles; ii++) {
190
191
192

        // multipoles

193
        int axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY;
194
195
196
        double charge, tholeD, dampingFactorD, polarityD;
        std::vector<double> dipolesD;
        std::vector<double> quadrupolesD;
197
        force.getMultipoleParameters(ii, charge, dipolesD, quadrupolesD, axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY,
198
                                     tholeD, dampingFactorD, polarityD);
199
200
201

        totalCharge                       += charge;
        axisTypes[ii]                      = axisType;
202
203
204
        multipoleAtomZs[ii]                = multipoleAtomZ;
        multipoleAtomXs[ii]                = multipoleAtomX;
        multipoleAtomYs[ii]                = multipoleAtomY;
205

peastman's avatar
peastman committed
206
207
208
209
        charges[ii]                        = charge;
        tholes[ii]                         = tholeD;
        dampingFactors[ii]                 = dampingFactorD;
        polarity[ii]                       = polarityD;
210

peastman's avatar
peastman committed
211
212
213
        dipoles[dipoleIndex++]             = dipolesD[0];
        dipoles[dipoleIndex++]             = dipolesD[1];
        dipoles[dipoleIndex++]             = dipolesD[2];
214
        
peastman's avatar
peastman committed
215
216
217
218
219
220
221
222
223
        quadrupoles[quadrupoleIndex++]     = quadrupolesD[0];
        quadrupoles[quadrupoleIndex++]     = quadrupolesD[1];
        quadrupoles[quadrupoleIndex++]     = quadrupolesD[2];
        quadrupoles[quadrupoleIndex++]     = quadrupolesD[3];
        quadrupoles[quadrupoleIndex++]     = quadrupolesD[4];
        quadrupoles[quadrupoleIndex++]     = quadrupolesD[5];
        quadrupoles[quadrupoleIndex++]     = quadrupolesD[6];
        quadrupoles[quadrupoleIndex++]     = quadrupolesD[7];
        quadrupoles[quadrupoleIndex++]     = quadrupolesD[8];
224
225
226
227

        // covalent info

        std::vector< std::vector<int> > covalentLists;
228
        force.getCovalentMaps(ii, covalentLists);
229
230
231
232
        multipoleAtomCovalentInfo[ii] = covalentLists;

    }

233
    polarizationType = force.getPolarizationType();
234
    if (polarizationType == AmoebaMultipoleForce::Mutual) {
235
236
        mutualInducedMaxIterations = force.getMutualInducedMaxIterations();
        mutualInducedTargetEpsilon = force.getMutualInducedTargetEpsilon();
237
238
    } else if (polarizationType == AmoebaMultipoleForce::Extrapolated) {
        extrapolationCoefficients = force.getExtrapolationCoefficients();
239
    }
240

241
242
243
    // PME

    nonbondedMethod  = force.getNonbondedMethod();
244
    if (nonbondedMethod == AmoebaMultipoleForce::PME) {
245
        usePme     = true;
peastman's avatar
peastman committed
246
247
        pmeGridDimension.resize(3);
        force.getPMEParameters(alphaEwald, pmeGridDimension[0], pmeGridDimension[1], pmeGridDimension[2]);
248
249
250
251
252
253
        cutoffDistance = force.getCutoffDistance();
        if (pmeGridDimension[0] == 0 || alphaEwald == 0.0) {
            NonbondedForce nb;
            nb.setEwaldErrorTolerance(force.getEwaldErrorTolerance());
            nb.setCutoffDistance(force.getCutoffDistance());
            int gridSizeX, gridSizeY, gridSizeZ;
254
            NonbondedForceImpl::calcPMEParameters(system, nb, alphaEwald, gridSizeX, gridSizeY, gridSizeZ, false);
255
256
257
258
259
260
261
262
            pmeGridDimension[0] = gridSizeX;
            pmeGridDimension[1] = gridSizeY;
            pmeGridDimension[2] = gridSizeZ;
        }    
    } else {
        usePme = false;
    }
    return;
263
264
}

265
AmoebaReferenceMultipoleForce* ReferenceCalcAmoebaMultipoleForceKernel::setupAmoebaReferenceMultipoleForce(ContextImpl& context)
266
{
267

268
269
270
271
272
    // amoebaReferenceMultipoleForce is set to AmoebaReferenceGeneralizedKirkwoodForce if AmoebaGeneralizedKirkwoodForce is present
    // amoebaReferenceMultipoleForce is set to AmoebaReferencePmeMultipoleForce if 'usePme' is set
    // amoebaReferenceMultipoleForce is set to AmoebaReferenceMultipoleForce otherwise

    // check if AmoebaGeneralizedKirkwoodForce is present 
273

274
275
276
277
278
279
280
281
282
    ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel* gkKernel = NULL;
    for (unsigned int ii = 0; ii < context.getForceImpls().size() && gkKernel == NULL; ii++) {
        AmoebaGeneralizedKirkwoodForceImpl* gkImpl = dynamic_cast<AmoebaGeneralizedKirkwoodForceImpl*>(context.getForceImpls()[ii]);
        if (gkImpl != NULL) {
            gkKernel = dynamic_cast<ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel*>(&gkImpl->getKernel().getImpl());
        }
    }    

    AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = NULL;
283
    if (gkKernel) {
284
285
286
287
288

        // amoebaReferenceGeneralizedKirkwoodForce is deleted in AmoebaReferenceGeneralizedKirkwoodMultipoleForce
        // destructor

        AmoebaReferenceGeneralizedKirkwoodForce* amoebaReferenceGeneralizedKirkwoodForce = new AmoebaReferenceGeneralizedKirkwoodForce();
289
290
291
292
293
294
295
296
        amoebaReferenceGeneralizedKirkwoodForce->setNumParticles(gkKernel->getNumParticles());
        amoebaReferenceGeneralizedKirkwoodForce->setSoluteDielectric(gkKernel->getSoluteDielectric());
        amoebaReferenceGeneralizedKirkwoodForce->setSolventDielectric(gkKernel->getSolventDielectric());
        amoebaReferenceGeneralizedKirkwoodForce->setDielectricOffset(gkKernel->getDielectricOffset());
        amoebaReferenceGeneralizedKirkwoodForce->setProbeRadius(gkKernel->getProbeRadius());
        amoebaReferenceGeneralizedKirkwoodForce->setSurfaceAreaFactor(gkKernel->getSurfaceAreaFactor());
        amoebaReferenceGeneralizedKirkwoodForce->setIncludeCavityTerm(gkKernel->getIncludeCavityTerm());
        amoebaReferenceGeneralizedKirkwoodForce->setDirectPolarization(gkKernel->getDirectPolarization());
297
298
299
300
301
302
303
304
305
306
        amoebaReferenceGeneralizedKirkwoodForce->setTanhRescaling(gkKernel->getTanhRescaling());
        double beta0, beta1, beta2;
        gkKernel->getTanhParameters(beta0, beta1, beta2);
        amoebaReferenceGeneralizedKirkwoodForce->setTanhParameters(beta0, beta1, beta2);
        amoebaReferenceGeneralizedKirkwoodForce->setDescreenOffset(gkKernel->getDescreenOffset());
        amoebaReferenceGeneralizedKirkwoodForce->setAtomicRadii(gkKernel->getAtomicRadii());
        amoebaReferenceGeneralizedKirkwoodForce->setScaleFactors(gkKernel->getScaleFactors());
        amoebaReferenceGeneralizedKirkwoodForce->setCharges(gkKernel->getCharges());
        amoebaReferenceGeneralizedKirkwoodForce->setDescreenRadii(gkKernel->getDescreenRadii());
        amoebaReferenceGeneralizedKirkwoodForce->setNeckFactors(gkKernel->getNeckFactors());
307
308

        // calculate Grycuk Born radii
peastman's avatar
peastman committed
309
        vector<Vec3>& posData   = extractPositions(context);
310
        amoebaReferenceGeneralizedKirkwoodForce->calculateGrycukBornRadii(posData);
311

312
        amoebaReferenceMultipoleForce = new AmoebaReferenceGeneralizedKirkwoodMultipoleForce(amoebaReferenceGeneralizedKirkwoodForce);
313

314
    } else if (usePme) {
315

316
317
318
319
        AmoebaReferencePmeMultipoleForce* amoebaReferencePmeMultipoleForce = new AmoebaReferencePmeMultipoleForce();
        amoebaReferencePmeMultipoleForce->setAlphaEwald(alphaEwald);
        amoebaReferencePmeMultipoleForce->setCutoffDistance(cutoffDistance);
        amoebaReferencePmeMultipoleForce->setPmeGridDimensions(pmeGridDimension);
peastman's avatar
peastman committed
320
        Vec3* boxVectors = extractBoxVectors(context);
321
        double minAllowedSize = 1.999999*cutoffDistance;
322
        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize) {
323
            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
324
325
326
        }
        amoebaReferencePmeMultipoleForce->setPeriodicBoxSize(boxVectors);
        amoebaReferenceMultipoleForce = static_cast<AmoebaReferenceMultipoleForce*>(amoebaReferencePmeMultipoleForce);
327

328
    } else {
329
         amoebaReferenceMultipoleForce = new AmoebaReferenceMultipoleForce(AmoebaReferenceMultipoleForce::NoCutoff);
330
331
    }

332
333
    // set polarization type

334
335
336
337
338
339
    if (polarizationType == AmoebaMultipoleForce::Mutual) {
        amoebaReferenceMultipoleForce->setPolarizationType(AmoebaReferenceMultipoleForce::Mutual);
        amoebaReferenceMultipoleForce->setMutualInducedDipoleTargetEpsilon(mutualInducedTargetEpsilon);
        amoebaReferenceMultipoleForce->setMaximumMutualInducedDipoleIterations(mutualInducedMaxIterations);
    } else if (polarizationType == AmoebaMultipoleForce::Direct) {
        amoebaReferenceMultipoleForce->setPolarizationType(AmoebaReferenceMultipoleForce::Direct);
340
341
342
    } else if (polarizationType == AmoebaMultipoleForce::Extrapolated) {
        amoebaReferenceMultipoleForce->setPolarizationType(AmoebaReferenceMultipoleForce::Extrapolated);
        amoebaReferenceMultipoleForce->setExtrapolationCoefficients(extrapolationCoefficients);
343
    } else {
344
        throw OpenMMException("Polarization type not recognzied.");
345
346
    }

347
348
349
350
351
352
    return amoebaReferenceMultipoleForce;

}

double ReferenceCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {

353
    AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
354

peastman's avatar
peastman committed
355
356
357
358
359
360
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy = amoebaReferenceMultipoleForce->calculateForceAndEnergy(posData, charges, dipoles, quadrupoles, tholes,
                                                                           dampingFactors, polarity, axisTypes, 
                                                                           multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
                                                                           multipoleAtomCovalentInfo, forceData);
361

362
    delete amoebaReferenceMultipoleForce;
363
364
365
366

    return static_cast<double>(energy);
}

367
368
369
370
371
372
void ReferenceCalcAmoebaMultipoleForceKernel::getInducedDipoles(ContextImpl& context, vector<Vec3>& outputDipoles) {
    int numParticles = context.getSystem().getNumParticles();
    outputDipoles.resize(numParticles);

    // Create an AmoebaReferenceMultipoleForce to do the calculation.
    
373
    AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
peastman's avatar
peastman committed
374
    vector<Vec3>& posData = extractPositions(context);
375
376
377
    
    // Retrieve the induced dipoles.
    
peastman's avatar
peastman committed
378
    vector<Vec3> inducedDipoles;
379
380
381
382
383
384
385
    amoebaReferenceMultipoleForce->calculateInducedDipoles(posData, charges, dipoles, quadrupoles, tholes,
            dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, multipoleAtomCovalentInfo, inducedDipoles);
    for (int i = 0; i < numParticles; i++)
        outputDipoles[i] = inducedDipoles[i];
    delete amoebaReferenceMultipoleForce;
}

386
387
388
389
390
391
392
void ReferenceCalcAmoebaMultipoleForceKernel::getLabFramePermanentDipoles(ContextImpl& context, vector<Vec3>& outputDipoles) {
    int numParticles = context.getSystem().getNumParticles();
    outputDipoles.resize(numParticles);

    // Create an AmoebaReferenceMultipoleForce to do the calculation.
    
    AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
peastman's avatar
peastman committed
393
    vector<Vec3>& posData = extractPositions(context);
394
395
396
    
    // Retrieve the permanent dipoles in the lab frame.
    
peastman's avatar
peastman committed
397
    vector<Vec3> labFramePermanentDipoles;
398
399
400
401
402
403
    amoebaReferenceMultipoleForce->calculateLabFramePermanentDipoles(posData, charges, dipoles, quadrupoles, tholes, 
            dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, multipoleAtomCovalentInfo, labFramePermanentDipoles);
    for (int i = 0; i < numParticles; i++)
        outputDipoles[i] = labFramePermanentDipoles[i];
    delete amoebaReferenceMultipoleForce;
}
404
405


406
407
408
409
410
411
412
void ReferenceCalcAmoebaMultipoleForceKernel::getTotalDipoles(ContextImpl& context, vector<Vec3>& outputDipoles) {
    int numParticles = context.getSystem().getNumParticles();
    outputDipoles.resize(numParticles);

    // Create an AmoebaReferenceMultipoleForce to do the calculation.
    
    AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
peastman's avatar
peastman committed
413
    vector<Vec3>& posData = extractPositions(context);
414
415
416
    
    // Retrieve the permanent dipoles in the lab frame.
    
peastman's avatar
peastman committed
417
    vector<Vec3> totalDipoles;
418
    amoebaReferenceMultipoleForce->calculateTotalDipoles(posData, charges, dipoles, quadrupoles, tholes,
419
            dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, multipoleAtomCovalentInfo, totalDipoles);
420
421
422
423
424
425
426

    for (int i = 0; i < numParticles; i++)
        outputDipoles[i] = totalDipoles[i];
    delete amoebaReferenceMultipoleForce;
}


427

428
void ReferenceCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
429
                                                                        std::vector< double >& outputElectrostaticPotential) {
430

431
    AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
peastman's avatar
peastman committed
432
433
434
    vector<Vec3>& posData                                     = extractPositions(context);
    vector<Vec3> grid(inputGrid.size());
    vector<double> potential(inputGrid.size());
435
    for (unsigned int ii = 0; ii < inputGrid.size(); ii++) {
436
437
        grid[ii] = inputGrid[ii];
    }
438
439
440
441
    amoebaReferenceMultipoleForce->calculateElectrostaticPotential(posData, charges, dipoles, quadrupoles, tholes,
                                                                   dampingFactors, polarity, axisTypes, 
                                                                   multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
                                                                   multipoleAtomCovalentInfo, grid, potential);
442

443
444
    outputElectrostaticPotential.resize(inputGrid.size());
    for (unsigned int ii = 0; ii < inputGrid.size(); ii++) {
445
446
447
448
        outputElectrostaticPotential[ii] = potential[ii];
    }

    delete amoebaReferenceMultipoleForce;
449
450
}

451
void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMoments) {
452
453
454

    // retrieve masses

455
    const System& system             = context.getSystem();
peastman's avatar
peastman committed
456
    vector<double> masses;
457
    for (int i = 0; i <  system.getNumParticles(); ++i) {
peastman's avatar
peastman committed
458
        masses.push_back(system.getParticleMass(i));
459
460
    }    

461
    AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
peastman's avatar
peastman committed
462
    vector<Vec3>& posData                                     = extractPositions(context);
463
464
465
466
    amoebaReferenceMultipoleForce->calculateAmoebaSystemMultipoleMoments(masses, posData, charges, dipoles, quadrupoles, tholes,
                                                                         dampingFactors, polarity, axisTypes, 
                                                                         multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
                                                                         multipoleAtomCovalentInfo, outputMultipoleMoments);
467
468

    delete amoebaReferenceMultipoleForce;
469
470
}

471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
void ReferenceCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaMultipoleForce& force) {
    if (numMultipoles != force.getNumMultipoles())
        throw OpenMMException("updateParametersInContext: The number of multipoles has changed");

    // Record the values.

    int dipoleIndex = 0;
    int quadrupoleIndex = 0;
    for (int i = 0; i < numMultipoles; ++i) {
        int axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY;
        double charge, tholeD, dampingFactorD, polarityD;
        std::vector<double> dipolesD;
        std::vector<double> quadrupolesD;
        force.getMultipoleParameters(i, charge, dipolesD, quadrupolesD, axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY, tholeD, dampingFactorD, polarityD);
        axisTypes[i] = axisType;
        multipoleAtomZs[i] = multipoleAtomZ;
        multipoleAtomXs[i] = multipoleAtomX;
        multipoleAtomYs[i] = multipoleAtomY;
peastman's avatar
peastman committed
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
        charges[i] = charge;
        tholes[i] = tholeD;
        dampingFactors[i] = dampingFactorD;
        polarity[i] = polarityD;
        dipoles[dipoleIndex++] = dipolesD[0];
        dipoles[dipoleIndex++] = dipolesD[1];
        dipoles[dipoleIndex++] = dipolesD[2];
        quadrupoles[quadrupoleIndex++] = quadrupolesD[0];
        quadrupoles[quadrupoleIndex++] = quadrupolesD[1];
        quadrupoles[quadrupoleIndex++] = quadrupolesD[2];
        quadrupoles[quadrupoleIndex++] = quadrupolesD[3];
        quadrupoles[quadrupoleIndex++] = quadrupolesD[4];
        quadrupoles[quadrupoleIndex++] = quadrupolesD[5];
        quadrupoles[quadrupoleIndex++] = quadrupolesD[6];
        quadrupoles[quadrupoleIndex++] = quadrupolesD[7];
        quadrupoles[quadrupoleIndex++] = quadrupolesD[8];
505
506
507
    }
}

508
509
510
511
512
513
514
515
516
void ReferenceCalcAmoebaMultipoleForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
    if (!usePme)
        throw OpenMMException("getPMEParametersInContext: This Context is not using PME");
    alpha = alphaEwald;
    nx = pmeGridDimension[0];
    ny = pmeGridDimension[1];
    nz = pmeGridDimension[2];
}

517
518
519
520
/* -------------------------------------------------------------------------- *
 *                       AmoebaGeneralizedKirkwood                            *
 * -------------------------------------------------------------------------- */

521
ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(const std::string& name, const Platform& platform, const System& system) :
522
523
524
525
526
527
           CalcAmoebaGeneralizedKirkwoodForceKernel(name, platform), system(system) {
}

ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::~ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel() {
}

528
int ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getNumParticles() const {
529
530
531
    return numParticles;
}

532
int ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getIncludeCavityTerm() const {
533
534
535
    return includeCavityTerm;
}

536
537
538
539
540
541
542
543
544
545
546
547
548
549
bool ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getTanhRescaling() const {
    return tanhRescaling;
}

void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getTanhParameters(double &b0, double &b1, double &b2) const {
    b0 = beta0;
    b1 = beta1;
    b2 = beta2;
}

double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDescreenOffset() const {
    return descreenOffset;
}

550
int ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDirectPolarization() const {
551
552
553
    return directPolarization;
}

peastman's avatar
peastman committed
554
double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSoluteDielectric() const {
555
556
557
    return soluteDielectric;
}

peastman's avatar
peastman committed
558
double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSolventDielectric() const {
559
560
561
    return solventDielectric;
}

peastman's avatar
peastman committed
562
double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDielectricOffset() const {
563
564
565
    return dielectricOffset;
}

peastman's avatar
peastman committed
566
double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getProbeRadius() const {
567
568
569
    return probeRadius;
}

peastman's avatar
peastman committed
570
double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getSurfaceAreaFactor() const {
571
572
573
    return surfaceAreaFactor;
}

574
575
576
577
578
579
580
581
582
583
const vector<double>& ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getAtomicRadii() const {
    return atomicRadii;
}

const vector<double>& ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getCharges() const {
    return charges;
}

const vector<double>& ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getScaleFactors() const {
    return scaleFactors;
584
585
}

586
587
const vector<double>& ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getDescreenRadii() const {
    return descreenRadii;
588
589
}

590
591
const vector<double>& ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::getNeckFactors() const {
    return neckFactors;
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
}

void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& system, const AmoebaGeneralizedKirkwoodForce& force) {

    // check that AmoebaMultipoleForce is present

    const AmoebaMultipoleForce* amoebaMultipoleForce = NULL;
    for (int ii = 0; ii < system.getNumForces() && amoebaMultipoleForce == NULL; ii++) {
        amoebaMultipoleForce = dynamic_cast<const AmoebaMultipoleForce*>(&system.getForce(ii));
    }

    if (amoebaMultipoleForce == NULL) {
        throw OpenMMException("AmoebaGeneralizedKirkwoodForce requires the System to also contain an AmoebaMultipoleForce.");
    }

607
    if (amoebaMultipoleForce->getNonbondedMethod() != AmoebaMultipoleForce::NoCutoff) {
608
609
610
611
612
        throw OpenMMException("AmoebaGeneralizedKirkwoodForce requires the AmoebaMultipoleForce use the NoCutoff nonbonded method.");
    }

    numParticles = system.getNumParticles();

613
    for (int ii = 0; ii < numParticles; ii++) {
614
615
        double particleCharge, particleRadius, scalingFactor, descreenRadius, neckFactor;
        force.getParticleParameters(ii, particleCharge, particleRadius, scalingFactor, descreenRadius, neckFactor);
peastman's avatar
peastman committed
616
617
618
        atomicRadii.push_back(particleRadius);
        scaleFactors.push_back(scalingFactor);
        charges.push_back(particleCharge);
619
620
        descreenRadii.push_back(descreenRadius);
        neckFactors.push_back(neckFactor);
621
622
623
624
625
        // Make sure the charge matches the one specified by the AmoebaMultipoleForce.

        double charge2, thole, damping, polarity;
        int axisType, atomX, atomY, atomZ;
        vector<double> dipole, quadrupole;
626
627
        amoebaMultipoleForce->getMultipoleParameters(ii, charge2, dipole, quadrupole, axisType, atomZ, atomX, atomY, thole, damping, polarity);
        if (particleCharge != charge2) {
628
629
630
            throw OpenMMException("AmoebaGeneralizedKirkwoodForce and AmoebaMultipoleForce must specify the same charge for every atom.");
        }

631
    }
632
    includeCavityTerm  = force.getIncludeCavityTerm();
633
634
635
636
637
638
639
    tanhRescaling      = force.getTanhRescaling();
    double b0, b1, b2;
    force.getTanhParameters(b0, b1, b2);
    beta0 = b0;
    beta1 = b1;
    beta2 = b2;
    descreenOffset     = force.getDescreenOffset();
peastman's avatar
peastman committed
640
641
    soluteDielectric   = force.getSoluteDielectric();
    solventDielectric  = force.getSolventDielectric();
642
643
644
    dielectricOffset   = force.getDielectricOffset();
    probeRadius        = force.getProbeRadius();
    surfaceAreaFactor  = force.getSurfaceAreaFactor();
645
646
647
648
649
650
651
    directPolarization = amoebaMultipoleForce->getPolarizationType() == AmoebaMultipoleForce::Direct ? 1 : 0;
}

double ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
    // handled in AmoebaReferenceGeneralizedKirkwoodMultipoleForce, a derived class of the class AmoebaReferenceMultipoleForce
    return 0.0;
}
Mark Friedrichs's avatar
Mark Friedrichs committed
652

653
654
655
656
657
658
659
void ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaGeneralizedKirkwoodForce& force) {
    if (numParticles != force.getNumParticles())
        throw OpenMMException("updateParametersInContext: The number of particles has changed");

    // Record the values.

    for (int i = 0; i < numParticles; ++i) {
660
661
        double particleCharge, particleRadius, scalingFactor, descreenRadius, neckFactor;
        force.getParticleParameters(i, particleCharge, particleRadius, scalingFactor, descreenRadius, neckFactor);
662
663
664
        atomicRadii[i] = particleRadius;
        scaleFactors[i] = scalingFactor;
        charges[i] = particleCharge;
665
666
        descreenRadii[i] = descreenRadius;
        neckFactors[i] = neckFactor;
667
668
669
    }
}

670
ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(const std::string& name, const Platform& platform, const System& system) :
Mark Friedrichs's avatar
Mark Friedrichs committed
671
       CalcAmoebaVdwForceKernel(name, platform), system(system) {
672
    useCutoff = 0;
673
674
    usePBC = 0;
    cutoff = 1.0e+10;
675
    neighborList = NULL;
Mark Friedrichs's avatar
Mark Friedrichs committed
676
677
678
}

ReferenceCalcAmoebaVdwForceKernel::~ReferenceCalcAmoebaVdwForceKernel() {
679
    if (neighborList)
680
        delete neighborList;
Mark Friedrichs's avatar
Mark Friedrichs committed
681
682
683
684
685
686
687
}

void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) {

    // per-particle parameters

    numParticles = system.getNumParticles();
688
689
    useCutoff              = (force.getNonbondedMethod() != AmoebaVdwForce::NoCutoff);
    usePBC                 = (force.getNonbondedMethod() == AmoebaVdwForce::CutoffPeriodic);
peastman's avatar
peastman committed
690
    cutoff                 = force.getCutoffDistance();
691
692
    neighborList           = useCutoff ? new NeighborList() : NULL;
    dispersionCoefficient  = force.getUseDispersionCorrection() ?  AmoebaVdwForceImpl::calcDispersionCorrection(system, force) : 0.0;
693
    vdwForce.initialize(force);
Mark Friedrichs's avatar
Mark Friedrichs committed
694
695
696
697
}

double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {

peastman's avatar
peastman committed
698
699
700
    vector<Vec3>& posData   = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    double energy;
701
    double lambda = context.getParameter(AmoebaVdwForce::Lambda());
702
    if (useCutoff) {
703
        computeNeighborListVoxelHash(*neighborList, numParticles, posData, vdwForce.getExclusions(), extractBoxVectors(context), usePBC, cutoff, 0.0);
704
        if (usePBC) {
peastman's avatar
peastman committed
705
            Vec3* boxVectors = extractBoxVectors(context);
706
            double minAllowedSize = 1.999999*cutoff;
707
            if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
708
                throw OpenMMException("The periodic box size has decreased to less than twice the cutoff.");
709
            vdwForce.setPeriodicBox(boxVectors);
710
            energy  = vdwForce.calculateForceAndEnergy(numParticles, lambda, posData, *neighborList, forceData);
711
            energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
712
        }
713
714
    }
    else
715
        energy = vdwForce.calculateForceAndEnergy(numParticles, lambda, posData, forceData);
Mark Friedrichs's avatar
Mark Friedrichs committed
716
717
718
    return static_cast<double>(energy);
}

719
720
721
void ReferenceCalcAmoebaVdwForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaVdwForce& force) {
    if (numParticles != force.getNumParticles())
        throw OpenMMException("updateParametersInContext: The number of particles has changed");
722
    vdwForce.initialize(force);
723
724
}

Mark Friedrichs's avatar
Mark Friedrichs committed
725
726
727
728
/* -------------------------------------------------------------------------- *
 *                           AmoebaWcaDispersion                              *
 * -------------------------------------------------------------------------- */

729
ReferenceCalcAmoebaWcaDispersionForceKernel::ReferenceCalcAmoebaWcaDispersionForceKernel(const std::string& name, const Platform& platform, const System& system) :
Mark Friedrichs's avatar
Mark Friedrichs committed
730
731
732
733
734
735
736
737
738
739
740
741
742
           CalcAmoebaWcaDispersionForceKernel(name, platform), system(system) {
}

ReferenceCalcAmoebaWcaDispersionForceKernel::~ReferenceCalcAmoebaWcaDispersionForceKernel() {
}

void ReferenceCalcAmoebaWcaDispersionForceKernel::initialize(const System& system, const AmoebaWcaDispersionForce& force) {

    // per-particle parameters

    numParticles = system.getNumParticles();
    radii.resize(numParticles);
    epsilons.resize(numParticles);
743
    for (int ii = 0; ii < numParticles; ii++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
744
745

        double radius, epsilon;
746
        force.getParticleParameters(ii, radius, epsilon);
Mark Friedrichs's avatar
Mark Friedrichs committed
747

peastman's avatar
peastman committed
748
749
        radii[ii] = radius;
        epsilons[ii] = epsilon;
Mark Friedrichs's avatar
Mark Friedrichs committed
750
751
    }   

peastman's avatar
peastman committed
752
    totalMaximumDispersionEnergy = AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(force);
Mark Friedrichs's avatar
Mark Friedrichs committed
753

peastman's avatar
peastman committed
754
755
756
757
758
759
760
761
    epso    = force.getEpso();
    epsh    = force.getEpsh();
    rmino   = force.getRmino();
    rminh   = force.getRminh();
    awater  = force.getAwater();
    shctd   = force.getShctd();
    dispoff = force.getDispoff();
    slevy   = force.getSlevy();
Mark Friedrichs's avatar
Mark Friedrichs committed
762
763
764
}

double ReferenceCalcAmoebaWcaDispersionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
peastman's avatar
peastman committed
765
766
    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
767
    AmoebaReferenceWcaDispersionForce amoebaReferenceWcaDispersionForce(epso, epsh, rmino, rminh, awater, shctd, dispoff, slevy);
peastman's avatar
peastman committed
768
    double energy = amoebaReferenceWcaDispersionForce.calculateForceAndEnergy(numParticles, posData, radii, epsilons, totalMaximumDispersionEnergy, forceData);
Mark Friedrichs's avatar
Mark Friedrichs committed
769
770
    return static_cast<double>(energy);
}
771
772
773
774
775
776
777
778
779
780

void ReferenceCalcAmoebaWcaDispersionForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaWcaDispersionForce& force) {
    if (numParticles != force.getNumParticles())
        throw OpenMMException("updateParametersInContext: The number of particles has changed");

    // Record the values.

    for (int i = 0; i < numParticles; ++i) {
        double radius, epsilon;
        force.getParticleParameters(i, radius, epsilon);
peastman's avatar
peastman committed
781
782
        radii[i] = radius;
        epsilons[i] = epsilon;
783
    }
peastman's avatar
peastman committed
784
    totalMaximumDispersionEnergy = AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(force);
785
}
peastman's avatar
peastman committed
786
787
788
789
790
791


/* -------------------------------------------------------------------------- *
 *                              HippoNonbonded                                *
 * -------------------------------------------------------------------------- */

792
ReferenceCalcHippoNonbondedForceKernel::ReferenceCalcHippoNonbondedForceKernel(const std::string& name, const Platform& platform, const System& system) :
peastman's avatar
peastman committed
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
         CalcHippoNonbondedForceKernel(name, platform), ixn(NULL) {
}

ReferenceCalcHippoNonbondedForceKernel::~ReferenceCalcHippoNonbondedForceKernel() {
    if (ixn != NULL)
        delete ixn;
}

void ReferenceCalcHippoNonbondedForceKernel::initialize(const System& system, const HippoNonbondedForce& force) {
    numParticles = force.getNumParticles();
    if (force.getNonbondedMethod() == HippoNonbondedForce::PME)
        ixn = new AmoebaReferencePmeHippoNonbondedForce(force, system);
    else
        ixn = new AmoebaReferenceHippoNonbondedForce(force);
}

void ReferenceCalcHippoNonbondedForceKernel::setupAmoebaReferenceHippoNonbondedForce(ContextImpl& context) {
    if (ixn->getNonbondedMethod() == HippoNonbondedForce::PME) {
        AmoebaReferencePmeHippoNonbondedForce* force = dynamic_cast<AmoebaReferencePmeHippoNonbondedForce*>(ixn);
        Vec3* boxVectors = extractBoxVectors(context);
        double minAllowedSize = 1.999999*force->getCutoffDistance();
        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
        force->setPeriodicBoxSize(boxVectors);
    }
}

double ReferenceCalcHippoNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {

    setupAmoebaReferenceHippoNonbondedForce(context);

    vector<Vec3>& posData = extractPositions(context);
    vector<Vec3>& forceData = extractForces(context);
    return ixn->calculateForceAndEnergy(posData, forceData);
}

void ReferenceCalcHippoNonbondedForceKernel::getInducedDipoles(ContextImpl& context, vector<Vec3>& outputDipoles) {
    outputDipoles.resize(numParticles);
    setupAmoebaReferenceHippoNonbondedForce(context);
    vector<Vec3>& posData = extractPositions(context);
    
    // Retrieve the induced dipoles.
    
    vector<Vec3> inducedDipoles;
    ixn->calculateInducedDipoles(posData, inducedDipoles);
    for (int i = 0; i < numParticles; i++)
        outputDipoles[i] = inducedDipoles[i];
}

void ReferenceCalcHippoNonbondedForceKernel::getLabFramePermanentDipoles(ContextImpl& context, vector<Vec3>& outputDipoles) {
    outputDipoles.resize(numParticles);
    setupAmoebaReferenceHippoNonbondedForce(context);
    vector<Vec3>& posData = extractPositions(context);
    
    // Retrieve the permanent dipoles in the lab frame.
    
    vector<Vec3> labFramePermanentDipoles;
    ixn->calculateLabFramePermanentDipoles(posData, labFramePermanentDipoles);
    for (int i = 0; i < numParticles; i++)
        outputDipoles[i] = labFramePermanentDipoles[i];
}

void ReferenceCalcHippoNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const HippoNonbondedForce& force) {
    if (numParticles != force.getNumParticles())
        throw OpenMMException("updateParametersInContext: The number of multipoles has changed");
    delete ixn;
    ixn = NULL;
    if (force.getNonbondedMethod() == HippoNonbondedForce::PME)
        ixn = new AmoebaReferencePmeHippoNonbondedForce(force, context.getSystem());
    else
        ixn = new AmoebaReferenceHippoNonbondedForce(force);
}

void ReferenceCalcHippoNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
    if (ixn->getNonbondedMethod() != HippoNonbondedForce::PME)
        throw OpenMMException("getPMEParametersInContext: This Context is not using PME");
    AmoebaReferencePmeHippoNonbondedForce* force = dynamic_cast<AmoebaReferencePmeHippoNonbondedForce*>(ixn);
    alpha = force->getAlphaEwald();
    vector<int> dim;
    force->getPmeGridDimensions(dim);
    nx = dim[0];
    ny = dim[1];
    nz = dim[2];
}

void ReferenceCalcHippoNonbondedForceKernel::getDPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
    if (ixn->getNonbondedMethod() != HippoNonbondedForce::PME)
        throw OpenMMException("getDPMEParametersInContext: This Context is not using PME");
    AmoebaReferencePmeHippoNonbondedForce* force = dynamic_cast<AmoebaReferencePmeHippoNonbondedForce*>(ixn);
    alpha = force->getDispersionAlphaEwald();
    vector<int> dim;
    force->getDispersionPmeGridDimensions(dim);
    nx = dim[0];
    ny = dim[1];
    nz = dim[2];
}