ReferenceLJCoulombIxn.cpp 26.2 KB
Newer Older
1

2
/* Portions copyright (c) 2006-2020 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), periodicExceptions(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
void ReferenceLJCoulombIxn::setPeriodicExceptions(bool periodic) {
    periodicExceptions = periodic;
}

174
175
176
177
178
179
180
/**---------------------------------------------------------------------------------------

   Calculate Ewald ixn

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

   --------------------------------------------------------------------------------------- */
189

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

peastman's avatar
peastman committed
195
    static const double epsilon     =  1.0;
196

197
    int kmax                            = (ewald ? std::max(numRx, std::max(numRy,numRz)) : 0);
peastman's avatar
peastman committed
198
199
200
201
    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;
202

peastman's avatar
peastman committed
203
204
205
    double totalSelfEwaldEnergy     = 0.0;
    double realSpaceEwaldEnergy     = 0.0;
    double recipEnergy              = 0.0;
peastman's avatar
peastman committed
206
    double recipDispersionEnergy    = 0.0;
peastman's avatar
peastman committed
207
208
    double totalRecipEnergy         = 0.0;
    double vdwEnergy                = 0.0;
209

210
211
212
213
214
215
216
217
218
    // 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
    // **************************************************************************************
219

220
    if (includeReciprocal) {
221
        double totalCharge = 0.0;
222
        for (int atomID = 0; atomID < numberOfAtoms; atomID++) {
223
224
225
            double selfEwaldEnergy = ONE_4PI_EPS0*atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI;
            totalCharge += atomParameters[atomID][QIndex];
            if (ljpme) {
226
227
228
                // Dispersion self term
                selfEwaldEnergy -= pow(alphaDispersionEwald, 6.0) * 64.0*pow(atomParameters[atomID][SigIndex], 6.0) * pow(atomParameters[atomID][EpsIndex], 2.0) / 12.0;
            }
229
            totalSelfEwaldEnergy -= selfEwaldEnergy;
230
        }
231
232
233
        // Correction for the neutralizing plasma.
        double volume = periodicBoxVectors[0][0] * periodicBoxVectors[1][1] * periodicBoxVectors[2][2];
        totalSelfEwaldEnergy -= totalCharge*totalCharge/(8*EPSILON0*volume*alphaEwald*alphaEwald);
234
235
    }

236
    if (totalEnergy) {
237
238
239
        *totalEnergy += totalSelfEwaldEnergy;
    }

240
241
242
    // **************************************************************************************
    // RECIPROCAL SPACE EWALD ENERGY AND FORCES
    // **************************************************************************************
243
    // PME
244

245
    if (pme && includeReciprocal) {
Peter Eastman's avatar
Peter Eastman committed
246
        ReferencePME pme(alphaEwald,numberOfAtoms,meshDim,5,1);
247

peastman's avatar
peastman committed
248
        vector<double> charges(numberOfAtoms);
249
250
        for (int i = 0; i < numberOfAtoms; i++)
            charges[i] = atomParameters[i][QIndex];
Peter Eastman's avatar
Peter Eastman committed
251
        pme.exec(atomCoordinates, forces, charges, periodicBoxVectors, recipEnergy);
252

253
254
        if (totalEnergy)
            *totalEnergy += recipEnergy;
255

256
257
        if (ljpme) {
            // Dispersion reciprocal space terms
Peter Eastman's avatar
Peter Eastman committed
258
            ReferencePME ljpme(alphaDispersionEwald,numberOfAtoms,dispersionMeshDim,5,1);
259

peastman's avatar
peastman committed
260
261
            std::vector<Vec3> dpmeforces(numberOfAtoms);
            for (int i = 0; i < numberOfAtoms; i++)
262
                charges[i] = 8.0*pow(atomParameters[i][SigIndex], 3.0) * atomParameters[i][EpsIndex];
Peter Eastman's avatar
Peter Eastman committed
263
            ljpme.exec_dpme(atomCoordinates, dpmeforces, charges, periodicBoxVectors, recipDispersionEnergy);
peastman's avatar
peastman committed
264
265
            for (int i = 0; i < numberOfAtoms; i++)
                forces[i] += dpmeforces[i];
266
267
268
269
            if (totalEnergy)
                *totalEnergy += recipDispersionEnergy;
        }
    }
270
271
    // Ewald method

272
    else if (ewald && includeReciprocal) {
273

274
        // setup reciprocal box
275

peastman's avatar
peastman committed
276
        double recipBoxSize[3] = { TWO_PI / periodicBoxVectors[0][0], TWO_PI / periodicBoxVectors[1][1], TWO_PI / periodicBoxVectors[2][2]};
277
278


279
        // setup K-vectors
280

281
282
283
284
#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);
285

286
287
        if (kmax < 1)
            throw OpenMMException("kmax for Ewald summation < 1");
288

289
290
291
        for (int i = 0; (i < numberOfAtoms); i++) {
            for (int m = 0; (m < 3); m++)
                EIR(0, i, m) = d_complex(1,0);
292

293
294
295
            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]));
296

297
298
299
300
            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);
        }
301

302
        // calculate reciprocal space energy and forces
303

304
305
        int lowry = 0;
        int lowrz = 1;
306

307
        for (int rx = 0; rx < numRx; rx++) {
308

peastman's avatar
peastman committed
309
            double kx = rx * recipBoxSize[0];
310

311
            for (int ry = lowry; ry < numRy; ry++) {
312

peastman's avatar
peastman committed
313
                double ky = ry * recipBoxSize[1];
314

315
316
317
318
                if (ry >= 0) {
                    for (int n = 0; n < numberOfAtoms; n++)
                        tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
                }
319

320
321
322
323
                else {
                    for (int n = 0; n < numberOfAtoms; n++)
                        tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
                }
324

325
                for (int rz = lowrz; rz < numRz; rz++) {
326

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

332
333
334
335
                    else {
                        for (int n = 0; n < numberOfAtoms; n++)
                            tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2)));
                    }
336

peastman's avatar
peastman committed
337
338
                    double cs = 0.0f;
                    double ss = 0.0f;
339

340
341
342
343
                    for (int n = 0; n < numberOfAtoms; n++) {
                        cs += tab_qxyz[n].real();
                        ss += tab_qxyz[n].imag();
                    }
344

peastman's avatar
peastman committed
345
346
347
                    double kz = rz * recipBoxSize[2];
                    double k2 = kx * kx + ky * ky + kz * kz;
                    double ak = exp(k2*factorEwald) / k2;
348

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

356
357
                    recipEnergy       = recipCoeff * ak * (cs * cs + ss * ss);
                    totalRecipEnergy += recipEnergy;
358

359
360
                    if (totalEnergy)
                        *totalEnergy += recipEnergy;
361

362
363
364
365
                    lowrz = 1 - numRz;
                }
                lowry = 1 - numRy;
            }
366
367
368
        }
    }

369
370
371
    // **************************************************************************************
    // SHORT-RANGE ENERGY AND FORCES
    // **************************************************************************************
372

373
374
    if (!includeDirect)
        return;
peastman's avatar
peastman committed
375
376
    double totalVdwEnergy            = 0.0f;
    double totalRealSpaceEwaldEnergy = 0.0f;
377

378

peastman's avatar
peastman committed
379
    for (auto& pair : *neighborList) {
380
381
382
        int ii = pair.first;
        int jj = pair.second;

peastman's avatar
peastman committed
383
        double deltaR[2][ReferenceForce::LastDeltaRIndex];
384
        ReferenceForce::getDeltaRPeriodic(atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0]);
peastman's avatar
peastman committed
385
386
387
        double r         = deltaR[0][ReferenceForce::RIndex];
        double inverseR  = 1.0/(deltaR[0][ReferenceForce::RIndex]);
        double switchValue = 1, switchDeriv = 0;
388
        if (useSwitch && r > switchingDistance) {
peastman's avatar
peastman committed
389
            double t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
390
391
392
            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
393
        double alphaR = alphaEwald * r;
394
395


peastman's avatar
peastman committed
396
397
        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);
398

peastman's avatar
peastman committed
399
400
401
402
403
404
405
        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;
406
407

        if (ljpme) {
peastman's avatar
peastman committed
408
409
410
411
412
413
414
            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];
415
416
417
418
            // 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
419
            double emult = c6i*c6j*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4));
420
            dEdR += 6.0*c6i*c6j*inverseR2*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4 + dar6/6.0));
421

peastman's avatar
peastman committed
422
423
            double inverseCut2 = 1.0/(cutoffDistance*cutoffDistance);
            double inverseCut6 = inverseCut2*inverseCut2*inverseCut2;
424
425
426
            sig2 = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
            sig2 *= sig2;
            sig6 = sig2*sig2*sig2;
427
            // The additive part of the potential shift
peastman's avatar
peastman committed
428
            double potentialshift = eps*(1.0-sig6*inverseCut6)*sig6*inverseCut6;
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
            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
445
            double force  = dEdR*deltaR[0][kk];
446
447
448
449
450
451
            forces[ii][kk]   += force;
            forces[jj][kk]   -= force;
        }

        // accumulate energies

peastman's avatar
peastman committed
452
        realSpaceEwaldEnergy        = ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR);
453
454
455

        totalVdwEnergy             += vdwEnergy;
        totalRealSpaceEwaldEnergy  += realSpaceEwaldEnergy;
456

457
    }
458

459
    if (totalEnergy)
460
461
        *totalEnergy += totalRealSpaceEwaldEnergy + totalVdwEnergy;

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

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

peastman's avatar
peastman committed
472
                double deltaR[2][ReferenceForce::LastDeltaRIndex];
473
474
475
476
                if (periodicExceptions)
                    ReferenceForce::getDeltaRPeriodic(atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0]);
                else
                    ReferenceForce::getDeltaR(atomCoordinates[jj], atomCoordinates[ii], deltaR[0]);
peastman's avatar
peastman committed
477
478
479
                double r         = deltaR[0][ReferenceForce::RIndex];
                double inverseR  = 1.0/(deltaR[0][ReferenceForce::RIndex]);
                double alphaR    = alphaEwald * r;
480
                if (erf(alphaR) > 1e-6) {
peastman's avatar
peastman committed
481
482
                    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);
483
484
485
486

                    // accumulate forces

                    for (int kk = 0; kk < 3; kk++) {
peastman's avatar
peastman committed
487
488
489
                        double force = dEdR*deltaR[0][kk];
                        forces[ii][kk] -= force;
                        forces[jj][kk] += force;
490
491
492
493
                    }

                    // accumulate energies

peastman's avatar
peastman committed
494
                    realSpaceEwaldEnergy = ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erf(alphaR);
495
496
                }
                else {
peastman's avatar
peastman committed
497
                    realSpaceEwaldEnergy = alphaEwald*TWO_OVER_SQRT_PI*ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex];
498
499
500
501
                }

                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
502
503
504
505
506
507
508
                    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];
509
                    realSpaceEwaldEnergy -= c6i*c6j*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4));
peastman's avatar
peastman committed
510
                    double dEdR = -6.0*c6i*c6j*inverseR2*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4 + dar6/6.0));
511
                    for (int kk = 0; kk < 3; kk++) {
peastman's avatar
peastman committed
512
513
514
                        double force = dEdR*deltaR[0][kk];
                        forces[ii][kk] -= force;
                        forces[jj][kk] += force;
515
516
517
518
                    }
                }

                totalExclusionEnergy += realSpaceEwaldEnergy;
Peter Eastman's avatar
Peter Eastman committed
519
            }
520
        }
Peter Eastman's avatar
Peter Eastman committed
521

522
    if (totalEnergy)
523
        *totalEnergy -= totalExclusionEnergy;
524
525
526
}


527
528
529
530
531
532
533
/**---------------------------------------------------------------------------------------

   Calculate LJ Coulomb pair ixn

   @param numberOfAtoms    number of atoms
   @param atomCoordinates  atom coordinates
   @param atomParameters   atom parameters                             atomParameters[atomIndex][paramterIndex]
534
535
   @param exclusions       atom exclusion indices
                           exclusions[atomIndex] contains the list of exclusions for that atom
536
537
   @param forces           force array (forces added)
   @param totalEnergy      total energy
538
539
   @param includeDirect      true if direct space interactions should be included
   @param includeReciprocal  true if reciprocal space interactions should be included
540
541

   --------------------------------------------------------------------------------------- */
542

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

547
    if (ewald || pme || ljpme) {
548
        calculateEwaldIxn(numberOfAtoms, atomCoordinates, atomParameters, exclusions, forces,
549
550
551
552
553
554
                          totalEnergy, includeDirect, includeReciprocal);
        return;
    }
    if (!includeDirect)
        return;
    if (cutoff) {
peastman's avatar
peastman committed
555
        for (auto& pair : *neighborList)
556
            calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, totalEnergy);
557
558
559
560
561
562
563
    }
    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())
564
                    calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, totalEnergy);
565
566
        }
    }
567
}
568

569
/**---------------------------------------------------------------------------------------
570

571
     Calculate LJ Coulomb pair ixn between two atoms
572

573
574
575
576
577
578
     @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
579

580
     --------------------------------------------------------------------------------------- */
581

peastman's avatar
peastman committed
582
void ReferenceLJCoulombIxn::calculateOneIxn(int ii, int jj, vector<Vec3>& atomCoordinates,
583
584
                                            vector<vector<double> >& atomParameters, vector<Vec3>& forces,
                                            double* totalEnergy) const {
peastman's avatar
peastman committed
585
    double deltaR[2][ReferenceForce::LastDeltaRIndex];
586
587
588
589

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

    if (periodic)
590
        ReferenceForce::getDeltaRPeriodic(atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0]);
591
    else
592
        ReferenceForce::getDeltaR(atomCoordinates[jj], atomCoordinates[ii], deltaR[0]);
593

peastman's avatar
peastman committed
594
595
596
    double r2        = deltaR[0][ReferenceForce::R2Index];
    double inverseR  = 1.0/(deltaR[0][ReferenceForce::RIndex]);
    double switchValue = 1, switchDeriv = 0;
597
    if (useSwitch) {
peastman's avatar
peastman committed
598
        double r = deltaR[0][ReferenceForce::RIndex];
599
        if (r > switchingDistance) {
peastman's avatar
peastman committed
600
            double t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
601
602
603
604
            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
605
606
607
608
    double sig = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
    double sig2 = inverseR*sig;
    sig2 *= sig2;
    double sig6 = sig2*sig2*sig2;
609

peastman's avatar
peastman committed
610
611
    double eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
    double dEdR = switchValue*eps*(12.0*sig6 - 6.0)*sig6;
612
    if (cutoff)
peastman's avatar
peastman committed
613
        dEdR += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
614
    else
peastman's avatar
peastman committed
615
        dEdR += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
616
    dEdR     *= inverseR*inverseR;
peastman's avatar
peastman committed
617
    double energy = eps*(sig6-1.0)*sig6;
618
619
620
621
622
    if (useSwitch) {
        dEdR -= energy*switchDeriv*inverseR;
        energy *= switchValue;
    }
    if (cutoff)
peastman's avatar
peastman committed
623
        energy += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf);
624
    else
peastman's avatar
peastman committed
625
        energy += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
626

627
628
    // accumulate forces

629
    for (int kk = 0; kk < 3; kk++) {
peastman's avatar
peastman committed
630
        double force  = dEdR*deltaR[0][kk];
631
632
        forces[ii][kk]   += force;
        forces[jj][kk]   -= force;
633
634
635
636
    }

    // accumulate energies

637
    if (totalEnergy)
638
639
        *totalEnergy += energy;
}
640