ReferenceLJCoulombIxn.cpp 25.9 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
            std::vector<Vec3> dpmeforces;
256
257
            for (int i = 0; i < numberOfAtoms; i++){
                charges[i] = 8.0*pow(atomParameters[i][SigIndex], 3.0) * atomParameters[i][EpsIndex];
peastman's avatar
peastman committed
258
                dpmeforces.push_back(Vec3());
259
260
261
262
263
264
265
266
267
268
269
270
            }
            pme_exec_dpme(pmedata,atomCoordinates,dpmeforces,charges,periodicBoxVectors,&recipDispersionEnergy);
            for (int i = 0; i < numberOfAtoms; i++){
                forces[i][0] -= 2.0*dpmeforces[i][0];
                forces[i][1] -= 2.0*dpmeforces[i][1];
                forces[i][2] -= 2.0*dpmeforces[i][2];
            }
            if (totalEnergy)
                *totalEnergy += recipDispersionEnergy;
            pme_destroy(pmedata);
        }
    }
271
272
    // Ewald method

273
    else if (ewald && includeReciprocal) {
274

275
        // setup reciprocal box
276

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


280
        // setup K-vectors
281

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

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

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

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

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

303
        // calculate reciprocal space energy and forces
304

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

379

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

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


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

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

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

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

        // accumulate energies

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

        totalVdwEnergy             += vdwEnergy;
        totalRealSpaceEwaldEnergy  += realSpaceEwaldEnergy;
457

458
    }
459

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

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

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

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

                    // accumulate forces

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

                    // accumulate energies

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

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

                totalExclusionEnergy += realSpaceEwaldEnergy;
Peter Eastman's avatar
Peter Eastman committed
517
            }
518
        }
Peter Eastman's avatar
Peter Eastman committed
519

520
    if (totalEnergy)
521
        *totalEnergy -= totalExclusionEnergy;
522
523
524
}


525
526
527
528
529
530
531
/**---------------------------------------------------------------------------------------

   Calculate LJ Coulomb pair ixn

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

   --------------------------------------------------------------------------------------- */
540

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

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

567
/**---------------------------------------------------------------------------------------
568

569
     Calculate LJ Coulomb pair ixn between two atoms
570

571
572
573
574
575
576
     @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
577

578
     --------------------------------------------------------------------------------------- */
579

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

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

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

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

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

625
626
    // accumulate forces

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

    // accumulate energies

635
    if (totalEnergy)
636
637
        *totalEnergy += energy;
}
638