ReferenceLJCoulombIxn.cpp 25.7 KB
Newer Older
1

2
/* Portions copyright (c) 2006-2018 Stanford University and Simbios.
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 * Contributors: Pande Group
 *
 * 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 <string.h>
#include <sstream>
27
#include <complex>
28
#include <algorithm>
29
#include <iostream>
30

31
#include "SimTKOpenMMUtilities.h"
32
33
#include "ReferenceLJCoulombIxn.h"
#include "ReferenceForce.h"
34
#include "ReferencePME.h"
peastman's avatar
peastman committed
35
#include "openmm/OpenMMException.h"
36

37
38
// In case we're using some primitive version of Visual Studio this will
// make sure that erf() and erfc() are defined.
39
#include "openmm/internal/MSVC_erfc.h"
40

41
using std::set;
42
using std::vector;
43
using namespace OpenMM;
44

45
46
47
48
49
50
/**---------------------------------------------------------------------------------------

   ReferenceLJCoulombIxn constructor

   --------------------------------------------------------------------------------------- */

51
ReferenceLJCoulombIxn::ReferenceLJCoulombIxn() : cutoff(false), useSwitch(false), periodic(false), ewald(false), pme(false), ljpme(false) {
52
53
54
55
56
57
58
59
}

/**---------------------------------------------------------------------------------------

   ReferenceLJCoulombIxn destructor

   --------------------------------------------------------------------------------------- */

60
ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn() {
61
62
}

63
/**---------------------------------------------------------------------------------------
64
65
66
67
68
69
70
71
72

     Set the force to use a cutoff.

     @param distance            the cutoff distance
     @param neighbors           the neighbor list to use
     @param solventDielectric   the dielectric constant of the bulk solvent

     --------------------------------------------------------------------------------------- */

peastman's avatar
peastman committed
73
void ReferenceLJCoulombIxn::setUseCutoff(double distance, const OpenMM::NeighborList& neighbors, double solventDielectric) {
74

75
76
77
    cutoff = true;
    cutoffDistance = distance;
    neighborList = &neighbors;
78
79
    krf = pow(cutoffDistance, -3.0)*(solventDielectric-1.0)/(2.0*solventDielectric+1.0);
    crf = (1.0/cutoffDistance)*(3.0*solventDielectric)/(2.0*solventDielectric+1.0);
80
}
81

82
83
84
85
86
87
88
89
/**---------------------------------------------------------------------------------------

   Set the force to use a switching function on the Lennard-Jones interaction.

   @param distance            the switching distance

   --------------------------------------------------------------------------------------- */

peastman's avatar
peastman committed
90
void ReferenceLJCoulombIxn::setUseSwitchingFunction(double distance) {
91
92
93
94
    useSwitch = true;
    switchingDistance = distance;
}

95
/**---------------------------------------------------------------------------------------
96
97
98
99
100

     Set the force to use periodic boundary conditions.  This requires that a cutoff has
     also been set, and the smallest side of the periodic box is at least twice the cutoff
     distance.

101
     @param vectors    the vectors defining the periodic box
102
103
104

     --------------------------------------------------------------------------------------- */

peastman's avatar
peastman committed
105
void ReferenceLJCoulombIxn::setPeriodic(OpenMM::Vec3* vectors) {
106
107

    assert(cutoff);
108
109
110
    assert(vectors[0][0] >= 2.0*cutoffDistance);
    assert(vectors[1][1] >= 2.0*cutoffDistance);
    assert(vectors[2][2] >= 2.0*cutoffDistance);
111
    periodic = true;
112
113
114
    periodicBoxVectors[0] = vectors[0];
    periodicBoxVectors[1] = vectors[1];
    periodicBoxVectors[2] = vectors[2];
115
}
116

117
/**---------------------------------------------------------------------------------------
118
119
120
121
122
123
124
125
126
127

     Set the force to use Ewald summation.

     @param alpha  the Ewald separation parameter
     @param kmaxx  the largest wave vector in the x direction
     @param kmaxy  the largest wave vector in the y direction
     @param kmaxz  the largest wave vector in the z direction

     --------------------------------------------------------------------------------------- */

peastman's avatar
peastman committed
128
void ReferenceLJCoulombIxn::setUseEwald(double alpha, int kmaxx, int kmaxy, int kmaxz) {
129
130
131
132
133
134
    alphaEwald = alpha;
    numRx = kmaxx;
    numRy = kmaxy;
    numRz = kmaxz;
    ewald = true;
}
135

136
/**---------------------------------------------------------------------------------------
137
138
139
140

     Set the force to use Particle-Mesh Ewald (PME) summation.

     @param alpha  the Ewald separation parameter
141
     @param gridSize the dimensions of the mesh
142
143
144

     --------------------------------------------------------------------------------------- */

peastman's avatar
peastman committed
145
void ReferenceLJCoulombIxn::setUsePME(double alpha, int meshSize[3]) {
146
147
148
149
150
151
152
153
154
    alphaEwald = alpha;
    meshDim[0] = meshSize[0];
    meshDim[1] = meshSize[1];
    meshDim[2] = meshSize[2];
    pme = true;
}

/**---------------------------------------------------------------------------------------

155
     Set the force to use Particle-Mesh Ewald (PME) summation for dispersion terms.
156

157
158
     @param alpha  the dispersion Ewald separation parameter
     @param gridSize the dimensions of the dispersion mesh
159
160
161

     --------------------------------------------------------------------------------------- */

peastman's avatar
peastman committed
162
void ReferenceLJCoulombIxn::setUseLJPME(double alpha, int meshSize[3]) {
163
164
165
166
    alphaDispersionEwald = alpha;
    dispersionMeshDim[0] = meshSize[0];
    dispersionMeshDim[1] = meshSize[1];
    dispersionMeshDim[2] = meshSize[2];
167
168
    ljpme = true;
}
169

170
171
172
173
174
175
176
/**---------------------------------------------------------------------------------------

   Calculate Ewald ixn

   @param numberOfAtoms    number of atoms
   @param atomCoordinates  atom coordinates
   @param atomParameters   atom parameters                             atomParameters[atomIndex][paramterIndex]
177
178
   @param exclusions       atom exclusion indices
                           exclusions[atomIndex] contains the list of exclusions for that atom
179
180
   @param forces           force array (forces added)
   @param totalEnergy      total energy
181
182
   @param includeDirect      true if direct space interactions should be included
   @param includeReciprocal  true if reciprocal space interactions should be included
183
184

   --------------------------------------------------------------------------------------- */
185

peastman's avatar
peastman committed
186
void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<Vec3>& atomCoordinates,
187
188
                                              vector<vector<double> >& atomParameters, vector<set<int> >& exclusions,
                                              vector<Vec3>& forces, double* totalEnergy, bool includeDirect, bool includeReciprocal) const {
peastman's avatar
peastman committed
189
    typedef std::complex<double> d_complex;
190

peastman's avatar
peastman committed
191
    static const double epsilon     =  1.0;
192

193
    int kmax                            = (ewald ? std::max(numRx, std::max(numRy,numRz)) : 0);
peastman's avatar
peastman committed
194
195
196
197
    double factorEwald              = -1 / (4*alphaEwald*alphaEwald);
    double SQRT_PI                  = sqrt(PI_M);
    double TWO_PI                   = 2.0 * PI_M;
    double recipCoeff               = ONE_4PI_EPS0*4*PI_M/(periodicBoxVectors[0][0] * periodicBoxVectors[1][1] * periodicBoxVectors[2][2]) /epsilon;
198

peastman's avatar
peastman committed
199
200
201
    double totalSelfEwaldEnergy     = 0.0;
    double realSpaceEwaldEnergy     = 0.0;
    double recipEnergy              = 0.0;
peastman's avatar
peastman committed
202
    double recipDispersionEnergy    = 0.0;
peastman's avatar
peastman committed
203
204
    double totalRecipEnergy         = 0.0;
    double vdwEnergy                = 0.0;
205

206
207
208
209
210
211
212
213
214
    // A couple of sanity checks for
    if(ljpme && useSwitch)
        throw OpenMMException("Switching cannot be used with LJPME");
    if(ljpme && !pme)
        throw OpenMMException("LJPME has been set, without PME being set");

    // **************************************************************************************
    // SELF ENERGY
    // **************************************************************************************
215

216
    if (includeReciprocal) {
217
        for (int atomID = 0; atomID < numberOfAtoms; atomID++) {
peastman's avatar
peastman committed
218
            double selfEwaldEnergy       = ONE_4PI_EPS0*atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI;
219
220
221
222
            if(ljpme) {
                // Dispersion self term
                selfEwaldEnergy -= pow(alphaDispersionEwald, 6.0) * 64.0*pow(atomParameters[atomID][SigIndex], 6.0) * pow(atomParameters[atomID][EpsIndex], 2.0) / 12.0;
            }
223
            totalSelfEwaldEnergy            -= selfEwaldEnergy;
224
        }
225
226
    }

227
    if (totalEnergy) {
228
229
230
        *totalEnergy += totalSelfEwaldEnergy;
    }

231
232
233
    // **************************************************************************************
    // RECIPROCAL SPACE EWALD ENERGY AND FORCES
    // **************************************************************************************
234
    // PME
235

236
237
    if (pme && includeReciprocal) {
        pme_t          pmedata; /* abstract handle for PME data */
238

239
        pme_init(&pmedata,alphaEwald,numberOfAtoms,meshDim,5,1);
240

peastman's avatar
peastman committed
241
        vector<double> charges(numberOfAtoms);
242
243
244
        for (int i = 0; i < numberOfAtoms; i++)
            charges[i] = atomParameters[i][QIndex];
        pme_exec(pmedata,atomCoordinates,forces,charges,periodicBoxVectors,&recipEnergy);
245

246
247
        if (totalEnergy)
            *totalEnergy += recipEnergy;
248

249
        pme_destroy(pmedata);
250

251
252
253
254
        if (ljpme) {
            // Dispersion reciprocal space terms
            pme_init(&pmedata,alphaDispersionEwald,numberOfAtoms,dispersionMeshDim,5,1);

peastman's avatar
peastman committed
255
256
            std::vector<Vec3> dpmeforces(numberOfAtoms);
            for (int i = 0; i < numberOfAtoms; i++)
257
258
                charges[i] = 8.0*pow(atomParameters[i][SigIndex], 3.0) * atomParameters[i][EpsIndex];
            pme_exec_dpme(pmedata,atomCoordinates,dpmeforces,charges,periodicBoxVectors,&recipDispersionEnergy);
peastman's avatar
peastman committed
259
260
            for (int i = 0; i < numberOfAtoms; i++)
                forces[i] += dpmeforces[i];
261
262
263
264
265
            if (totalEnergy)
                *totalEnergy += recipDispersionEnergy;
            pme_destroy(pmedata);
        }
    }
266
267
    // Ewald method

268
    else if (ewald && includeReciprocal) {
269

270
        // setup reciprocal box
271

peastman's avatar
peastman committed
272
        double recipBoxSize[3] = { TWO_PI / periodicBoxVectors[0][0], TWO_PI / periodicBoxVectors[1][1], TWO_PI / periodicBoxVectors[2][2]};
273
274


275
        // setup K-vectors
276

277
278
279
280
#define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
        vector<d_complex> eir(kmax*numberOfAtoms*3);
        vector<d_complex> tab_xy(numberOfAtoms);
        vector<d_complex> tab_qxyz(numberOfAtoms);
281

282
283
        if (kmax < 1)
            throw OpenMMException("kmax for Ewald summation < 1");
284

285
286
287
        for (int i = 0; (i < numberOfAtoms); i++) {
            for (int m = 0; (m < 3); m++)
                EIR(0, i, m) = d_complex(1,0);
288

289
290
291
            for (int m=0; (m<3); m++)
                EIR(1, i, m) = d_complex(cos(atomCoordinates[i][m]*recipBoxSize[m]),
                                         sin(atomCoordinates[i][m]*recipBoxSize[m]));
292

293
294
295
296
            for (int j=2; (j<kmax); j++)
                for (int m=0; (m<3); m++)
                    EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
        }
297

298
        // calculate reciprocal space energy and forces
299

300
301
        int lowry = 0;
        int lowrz = 1;
302

303
        for (int rx = 0; rx < numRx; rx++) {
304

peastman's avatar
peastman committed
305
            double kx = rx * recipBoxSize[0];
306

307
            for (int ry = lowry; ry < numRy; ry++) {
308

peastman's avatar
peastman committed
309
                double ky = ry * recipBoxSize[1];
310

311
312
313
314
                if (ry >= 0) {
                    for (int n = 0; n < numberOfAtoms; n++)
                        tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
                }
315

316
317
318
319
                else {
                    for (int n = 0; n < numberOfAtoms; n++)
                        tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
                }
320

321
                for (int rz = lowrz; rz < numRz; rz++) {
322

323
324
325
326
                    if (rz >= 0) {
                        for (int n = 0; n < numberOfAtoms; n++)
                            tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2));
                    }
327

328
329
330
331
                    else {
                        for (int n = 0; n < numberOfAtoms; n++)
                            tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2)));
                    }
332

peastman's avatar
peastman committed
333
334
                    double cs = 0.0f;
                    double ss = 0.0f;
335

336
337
338
339
                    for (int n = 0; n < numberOfAtoms; n++) {
                        cs += tab_qxyz[n].real();
                        ss += tab_qxyz[n].imag();
                    }
340

peastman's avatar
peastman committed
341
342
343
                    double kz = rz * recipBoxSize[2];
                    double k2 = kx * kx + ky * ky + kz * kz;
                    double ak = exp(k2*factorEwald) / k2;
344

345
                    for (int n = 0; n < numberOfAtoms; n++) {
peastman's avatar
peastman committed
346
                        double force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real());
347
348
349
350
                        forces[n][0] += 2 * recipCoeff * force * kx ;
                        forces[n][1] += 2 * recipCoeff * force * ky ;
                        forces[n][2] += 2 * recipCoeff * force * kz ;
                    }
351

352
353
                    recipEnergy       = recipCoeff * ak * (cs * cs + ss * ss);
                    totalRecipEnergy += recipEnergy;
354

355
356
                    if (totalEnergy)
                        *totalEnergy += recipEnergy;
357

358
359
360
361
                    lowrz = 1 - numRz;
                }
                lowry = 1 - numRy;
            }
362
363
364
        }
    }

365
366
367
    // **************************************************************************************
    // SHORT-RANGE ENERGY AND FORCES
    // **************************************************************************************
368

369
370
    if (!includeDirect)
        return;
peastman's avatar
peastman committed
371
372
    double totalVdwEnergy            = 0.0f;
    double totalRealSpaceEwaldEnergy = 0.0f;
373

374

peastman's avatar
peastman committed
375
    for (auto& pair : *neighborList) {
376
377
378
        int ii = pair.first;
        int jj = pair.second;

peastman's avatar
peastman committed
379
        double deltaR[2][ReferenceForce::LastDeltaRIndex];
380
        ReferenceForce::getDeltaRPeriodic(atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0]);
peastman's avatar
peastman committed
381
382
383
        double r         = deltaR[0][ReferenceForce::RIndex];
        double inverseR  = 1.0/(deltaR[0][ReferenceForce::RIndex]);
        double switchValue = 1, switchDeriv = 0;
384
        if (useSwitch && r > switchingDistance) {
peastman's avatar
peastman committed
385
            double t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
386
387
388
            switchValue = 1+t*t*t*(-10+t*(15-t*6));
            switchDeriv = t*t*(-30+t*(60-t*30))/(cutoffDistance-switchingDistance);
        }
peastman's avatar
peastman committed
389
        double alphaR = alphaEwald * r;
390
391


peastman's avatar
peastman committed
392
393
        double dEdR = ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
        dEdR = dEdR * (erfc(alphaR) + 2 * alphaR * exp (- alphaR * alphaR) / SQRT_PI);
394

peastman's avatar
peastman committed
395
396
397
398
399
400
401
        double sig = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
        double sig2 = inverseR*sig;
        sig2 *= sig2;
        double sig6 = sig2*sig2*sig2;
        double eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
        dEdR += switchValue*eps*(12.0*sig6 - 6.0)*sig6*inverseR*inverseR;
        vdwEnergy = eps*(sig6-1.0)*sig6;
402
403

        if (ljpme) {
peastman's avatar
peastman committed
404
405
406
407
408
409
410
            double dalphaR   = alphaDispersionEwald * r;
            double dar2 = dalphaR*dalphaR;
            double dar4 = dar2*dar2;
            double dar6 = dar4*dar2;
            double inverseR2 = inverseR*inverseR;
            double c6i = 8.0*pow(atomParameters[ii][SigIndex], 3.0) * atomParameters[ii][EpsIndex];
            double c6j = 8.0*pow(atomParameters[jj][SigIndex], 3.0) * atomParameters[jj][EpsIndex];
411
412
413
414
            // For the energies and forces, we first add the regular Lorentz−Berthelot terms.  The C12 term is treated as usual
            // but we then subtract out (remembering that the C6 term is negative) the multiplicative C6 term that has been
            // computed in real space.  Finally, we add a potential shift term to account for the difference between the LB
            // and multiplicative functional forms at the cutoff.
peastman's avatar
peastman committed
415
            double emult = c6i*c6j*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4));
416
            dEdR += 6.0*c6i*c6j*inverseR2*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4 + dar6/6.0));
417

peastman's avatar
peastman committed
418
419
            double inverseCut2 = 1.0/(cutoffDistance*cutoffDistance);
            double inverseCut6 = inverseCut2*inverseCut2*inverseCut2;
420
421
422
            sig2 = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
            sig2 *= sig2;
            sig6 = sig2*sig2*sig2;
423
            // The additive part of the potential shift
peastman's avatar
peastman committed
424
            double potentialshift = eps*(1.0-sig6*inverseCut6)*sig6*inverseCut6;
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
            dalphaR   = alphaDispersionEwald * cutoffDistance;
            dar2 = dalphaR*dalphaR;
            dar4 = dar2*dar2;
            // The multiplicative part of the potential shift
            potentialshift -= c6i*c6j*inverseCut6*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4));
            vdwEnergy += emult + potentialshift;
        }

        if (useSwitch) {
            dEdR -= vdwEnergy*switchDeriv*inverseR;
            vdwEnergy *= switchValue;
        }

        // accumulate forces

        for (int kk = 0; kk < 3; kk++) {
peastman's avatar
peastman committed
441
            double force  = dEdR*deltaR[0][kk];
442
443
444
445
446
447
            forces[ii][kk]   += force;
            forces[jj][kk]   -= force;
        }

        // accumulate energies

peastman's avatar
peastman committed
448
        realSpaceEwaldEnergy        = ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR);
449
450
451

        totalVdwEnergy             += vdwEnergy;
        totalRealSpaceEwaldEnergy  += realSpaceEwaldEnergy;
452

453
    }
454

455
    if (totalEnergy)
456
457
        *totalEnergy += totalRealSpaceEwaldEnergy + totalVdwEnergy;

Peter Eastman's avatar
Peter Eastman committed
458
459
    // Now subtract off the exclusions, since they were implicitly included in the reciprocal space sum.

peastman's avatar
peastman committed
460
    double totalExclusionEnergy = 0.0f;
461
    const double TWO_OVER_SQRT_PI = 2/sqrt(PI_M);
Peter Eastman's avatar
Peter Eastman committed
462
    for (int i = 0; i < numberOfAtoms; i++)
peastman's avatar
peastman committed
463
464
        for (int exclusion : exclusions[i]) {
            if (exclusion > i) {
465
                int ii = i;
peastman's avatar
peastman committed
466
                int jj = exclusion;
467

peastman's avatar
peastman committed
468
                double deltaR[2][ReferenceForce::LastDeltaRIndex];
469
                ReferenceForce::getDeltaR(atomCoordinates[jj], atomCoordinates[ii], deltaR[0]);
peastman's avatar
peastman committed
470
471
472
                double r         = deltaR[0][ReferenceForce::RIndex];
                double inverseR  = 1.0/(deltaR[0][ReferenceForce::RIndex]);
                double alphaR    = alphaEwald * r;
473
                if (erf(alphaR) > 1e-6) {
peastman's avatar
peastman committed
474
475
                    double dEdR = ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
                    dEdR = dEdR * (erf(alphaR) - 2 * alphaR * exp (- alphaR * alphaR) / SQRT_PI);
476
477
478
479

                    // accumulate forces

                    for (int kk = 0; kk < 3; kk++) {
peastman's avatar
peastman committed
480
481
482
                        double force = dEdR*deltaR[0][kk];
                        forces[ii][kk] -= force;
                        forces[jj][kk] += force;
483
484
485
486
                    }

                    // accumulate energies

peastman's avatar
peastman committed
487
                    realSpaceEwaldEnergy = ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erf(alphaR);
488
489
                }
                else {
peastman's avatar
peastman committed
490
                    realSpaceEwaldEnergy = alphaEwald*TWO_OVER_SQRT_PI*ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex];
491
492
493
494
                }

                if(ljpme){
                    // Dispersion terms.  Here we just back out the reciprocal space terms, and don't add any extra real space terms.
peastman's avatar
peastman committed
495
496
497
498
499
500
501
                    double dalphaR   = alphaDispersionEwald * r;
                    double inverseR2 = inverseR*inverseR;
                    double dar2 = dalphaR*dalphaR;
                    double dar4 = dar2*dar2;
                    double dar6 = dar4*dar2;
                    double c6i = 8.0*pow(atomParameters[ii][SigIndex], 3.0) * atomParameters[ii][EpsIndex];
                    double c6j = 8.0*pow(atomParameters[jj][SigIndex], 3.0) * atomParameters[jj][EpsIndex];
502
                    realSpaceEwaldEnergy -= c6i*c6j*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4));
peastman's avatar
peastman committed
503
                    double dEdR = -6.0*c6i*c6j*inverseR2*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4 + dar6/6.0));
504
                    for (int kk = 0; kk < 3; kk++) {
peastman's avatar
peastman committed
505
506
507
                        double force = dEdR*deltaR[0][kk];
                        forces[ii][kk] -= force;
                        forces[jj][kk] += force;
508
509
510
511
                    }
                }

                totalExclusionEnergy += realSpaceEwaldEnergy;
Peter Eastman's avatar
Peter Eastman committed
512
            }
513
        }
Peter Eastman's avatar
Peter Eastman committed
514

515
    if (totalEnergy)
516
        *totalEnergy -= totalExclusionEnergy;
517
518
519
}


520
521
522
523
524
525
526
/**---------------------------------------------------------------------------------------

   Calculate LJ Coulomb pair ixn

   @param numberOfAtoms    number of atoms
   @param atomCoordinates  atom coordinates
   @param atomParameters   atom parameters                             atomParameters[atomIndex][paramterIndex]
527
528
   @param exclusions       atom exclusion indices
                           exclusions[atomIndex] contains the list of exclusions for that atom
529
530
   @param forces           force array (forces added)
   @param totalEnergy      total energy
531
532
   @param includeDirect      true if direct space interactions should be included
   @param includeReciprocal  true if reciprocal space interactions should be included
533
534

   --------------------------------------------------------------------------------------- */
535

peastman's avatar
peastman committed
536
void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<Vec3>& atomCoordinates,
537
538
                                             vector<vector<double> >& atomParameters, vector<set<int> >& exclusions,
                                             vector<Vec3>& forces, double* totalEnergy, bool includeDirect, bool includeReciprocal) const {
539

540
    if (ewald || pme || ljpme) {
541
        calculateEwaldIxn(numberOfAtoms, atomCoordinates, atomParameters, exclusions, forces,
542
543
544
545
546
547
                          totalEnergy, includeDirect, includeReciprocal);
        return;
    }
    if (!includeDirect)
        return;
    if (cutoff) {
peastman's avatar
peastman committed
548
        for (auto& pair : *neighborList)
549
            calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, totalEnergy);
550
551
552
553
554
555
556
    }
    else {
        for (int ii = 0; ii < numberOfAtoms; ii++) {
            // loop over atom pairs

            for (int jj = ii+1; jj < numberOfAtoms; jj++)
                if (exclusions[jj].find(ii) == exclusions[jj].end())
557
                    calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, totalEnergy);
558
559
        }
    }
560
}
561

562
/**---------------------------------------------------------------------------------------
563

564
     Calculate LJ Coulomb pair ixn between two atoms
565

566
567
568
569
570
571
     @param ii               the index of the first atom
     @param jj               the index of the second atom
     @param atomCoordinates  atom coordinates
     @param atomParameters   atom parameters (charges, c6, c12, ...)     atomParameters[atomIndex][paramterIndex]
     @param forces           force array (forces added)
     @param totalEnergy      total energy
572

573
     --------------------------------------------------------------------------------------- */
574

peastman's avatar
peastman committed
575
void ReferenceLJCoulombIxn::calculateOneIxn(int ii, int jj, vector<Vec3>& atomCoordinates,
576
577
                                            vector<vector<double> >& atomParameters, vector<Vec3>& forces,
                                            double* totalEnergy) const {
peastman's avatar
peastman committed
578
    double deltaR[2][ReferenceForce::LastDeltaRIndex];
579
580
581
582

    // get deltaR, R2, and R between 2 atoms

    if (periodic)
583
        ReferenceForce::getDeltaRPeriodic(atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0]);
584
    else
585
        ReferenceForce::getDeltaR(atomCoordinates[jj], atomCoordinates[ii], deltaR[0]);
586

peastman's avatar
peastman committed
587
588
589
    double r2        = deltaR[0][ReferenceForce::R2Index];
    double inverseR  = 1.0/(deltaR[0][ReferenceForce::RIndex]);
    double switchValue = 1, switchDeriv = 0;
590
    if (useSwitch) {
peastman's avatar
peastman committed
591
        double r = deltaR[0][ReferenceForce::RIndex];
592
        if (r > switchingDistance) {
peastman's avatar
peastman committed
593
            double t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
594
595
596
597
            switchValue = 1+t*t*t*(-10+t*(15-t*6));
            switchDeriv = t*t*(-30+t*(60-t*30))/(cutoffDistance-switchingDistance);
        }
    }
peastman's avatar
peastman committed
598
599
600
601
    double sig = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
    double sig2 = inverseR*sig;
    sig2 *= sig2;
    double sig6 = sig2*sig2*sig2;
602

peastman's avatar
peastman committed
603
604
    double eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
    double dEdR = switchValue*eps*(12.0*sig6 - 6.0)*sig6;
605
    if (cutoff)
peastman's avatar
peastman committed
606
        dEdR += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
607
    else
peastman's avatar
peastman committed
608
        dEdR += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
609
    dEdR     *= inverseR*inverseR;
peastman's avatar
peastman committed
610
    double energy = eps*(sig6-1.0)*sig6;
611
612
613
614
615
    if (useSwitch) {
        dEdR -= energy*switchDeriv*inverseR;
        energy *= switchValue;
    }
    if (cutoff)
peastman's avatar
peastman committed
616
        energy += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf);
617
    else
peastman's avatar
peastman committed
618
        energy += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
619

620
621
    // accumulate forces

622
    for (int kk = 0; kk < 3; kk++) {
peastman's avatar
peastman committed
623
        double force  = dEdR*deltaR[0][kk];
624
625
        forces[ii][kk]   += force;
        forces[jj][kk]   -= force;
626
627
628
629
    }

    // accumulate energies

630
    if (totalEnergy)
631
632
        *totalEnergy += energy;
}
633