ReferenceLJCoulombIxn.cpp 21.7 KB
Newer Older
1

2
/* Portions copyright (c) 2006-2013 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

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

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

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

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

   ReferenceLJCoulombIxn constructor

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

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

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

   ReferenceLJCoulombIxn destructor

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

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

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

     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
72
  void ReferenceLJCoulombIxn::setUseCutoff(double distance, const OpenMM::NeighborList& neighbors, double solventDielectric) {
73

74
75
76
    cutoff = true;
    cutoffDistance = distance;
    neighborList = &neighbors;
77
78
    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);
79
80
  }

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

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

   @param distance            the switching distance

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

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

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

     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.

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

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

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

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

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

     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
127
  void ReferenceLJCoulombIxn::setUseEwald(double alpha, int kmaxx, int kmaxy, int kmaxz) {
128
129
130
131
132
133
134
      alphaEwald = alpha;
      numRx = kmaxx;
      numRy = kmaxy;
      numRz = kmaxz;
      ewald = true;
  }

135
136
137
138
139
  /**---------------------------------------------------------------------------------------

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

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

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

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

152
153
154
155
156
157
158
/**---------------------------------------------------------------------------------------

   Calculate Ewald ixn

   @param numberOfAtoms    number of atoms
   @param atomCoordinates  atom coordinates
   @param atomParameters   atom parameters                             atomParameters[atomIndex][paramterIndex]
159
160
   @param exclusions       atom exclusion indices
                           exclusions[atomIndex] contains the list of exclusions for that atom
161
162
163
164
   @param fixedParameters  non atom parameters (not currently used)
   @param forces           force array (forces added)
   @param energyByAtom     atom energy
   @param totalEnergy      total energy
165
166
   @param includeDirect      true if direct space interactions should be included
   @param includeReciprocal  true if reciprocal space interactions should be included
167
168

   --------------------------------------------------------------------------------------- */
169

peastman's avatar
peastman committed
170
171
172
173
174
void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<Vec3>& atomCoordinates,
                                             double** atomParameters, vector<set<int> >& exclusions,
                                             double* fixedParameters, vector<Vec3>& forces,
                                             double* energyByAtom, double* totalEnergy, bool includeDirect, bool includeReciprocal) const {
    typedef std::complex<double> d_complex;
175

peastman's avatar
peastman committed
176
    static const double epsilon     =  1.0;
177

178
    int kmax                            = (ewald ? std::max(numRx, std::max(numRy,numRz)) : 0);
peastman's avatar
peastman committed
179
180
181
182
    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;
183

peastman's avatar
peastman committed
184
185
186
187
188
    double totalSelfEwaldEnergy     = 0.0;
    double realSpaceEwaldEnergy     = 0.0;
    double recipEnergy              = 0.0;
    double totalRecipEnergy         = 0.0;
    double vdwEnergy                = 0.0;
189
190
191
192
193

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

194
    if (includeReciprocal) {
195
        for (int atomID = 0; atomID < numberOfAtoms; atomID++) {
peastman's avatar
peastman committed
196
            double selfEwaldEnergy       = ONE_4PI_EPS0*atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI;
197
            totalSelfEwaldEnergy            -= selfEwaldEnergy;
198
            if (energyByAtom) {
199
200
                energyByAtom[atomID]        -= selfEwaldEnergy;
            }
201
        }
202
203
    }

204
    if (totalEnergy) {
205
206
207
        *totalEnergy += totalSelfEwaldEnergy;
    }

208
209
210
// **************************************************************************************
// RECIPROCAL SPACE EWALD ENERGY AND FORCES
// **************************************************************************************
211
    // PME
212

213
  if (pme && includeReciprocal) {
214
    pme_t          pmedata; /* abstract handle for PME data */
215

216
    pme_init(&pmedata,alphaEwald,numberOfAtoms,meshDim,5,1);
217

peastman's avatar
peastman committed
218
    vector<double> charges(numberOfAtoms);
peastman's avatar
peastman committed
219
220
    for (int i = 0; i < numberOfAtoms; i++)
        charges[i] = atomParameters[i][QIndex];
221
    pme_exec(pmedata,atomCoordinates,forces,charges,periodicBoxVectors,&recipEnergy);
222

223
    if (totalEnergy)
224
225
       *totalEnergy += recipEnergy;

226
227
    if (energyByAtom)
        for (int n = 0; n < numberOfAtoms; n++)
228
229
            energyByAtom[n] += recipEnergy;

230
        pme_destroy(pmedata);
231
232
233
234
  }

    // Ewald method

235
  else if (ewald && includeReciprocal) {
236
237

    // setup reciprocal box
238

peastman's avatar
peastman committed
239
         double recipBoxSize[3] = { TWO_PI / periodicBoxVectors[0][0], TWO_PI / periodicBoxVectors[1][1], TWO_PI / periodicBoxVectors[2][2]};
240
241


242
    // setup K-vectors
243

244
  #define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
245
246
247
  vector<d_complex> eir(kmax*numberOfAtoms*3);
  vector<d_complex> tab_xy(numberOfAtoms);
  vector<d_complex> tab_qxyz(numberOfAtoms);
248

peastman's avatar
peastman committed
249
250
  if (kmax < 1)
      throw OpenMMException("kmax for Ewald summation < 1");
251

252
253
  for (int i = 0; (i < numberOfAtoms); i++) {
    for (int m = 0; (m < 3); m++)
254
      EIR(0, i, m) = d_complex(1,0);
255

256
    for (int m=0; (m<3); m++)
257
258
      EIR(1, i, m) = d_complex(cos(atomCoordinates[i][m]*recipBoxSize[m]),
                               sin(atomCoordinates[i][m]*recipBoxSize[m]));
259

260
261
    for (int j=2; (j<kmax); j++)
      for (int m=0; (m<3); m++)
262
        EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
263
264
  }

265
266
    // calculate reciprocal space energy and forces

267
268
269
    int lowry = 0;
    int lowrz = 1;

270
    for (int rx = 0; rx < numRx; rx++) {
271

peastman's avatar
peastman committed
272
      double kx = rx * recipBoxSize[0];
273

274
      for (int ry = lowry; ry < numRy; ry++) {
275

peastman's avatar
peastman committed
276
        double ky = ry * recipBoxSize[1];
277

278
279
        if (ry >= 0) {
          for (int n = 0; n < numberOfAtoms; n++)
280
            tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
281
282
283
        }

        else {
284
          for (int n = 0; n < numberOfAtoms; n++)
285
            tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
286
287
288
289
        }

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

290
291
          if (rz >= 0) {
           for (int n = 0; n < numberOfAtoms; n++)
292
             tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2));
293
294
295
          }

          else {
296
            for (int n = 0; n < numberOfAtoms; n++)
297
              tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2)));
298
299
          }

peastman's avatar
peastman committed
300
301
          double cs = 0.0f;
          double ss = 0.0f;
302

303
          for (int n = 0; n < numberOfAtoms; n++) {
304
305
306
307
            cs += tab_qxyz[n].real();
            ss += tab_qxyz[n].imag();
          }

peastman's avatar
peastman committed
308
309
310
          double kz = rz * recipBoxSize[2];
          double k2 = kx * kx + ky * ky + kz * kz;
          double ak = exp(k2*factorEwald) / k2;
311

312
          for (int n = 0; n < numberOfAtoms; n++) {
peastman's avatar
peastman committed
313
            double force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real());
314
315
316
            forces[n][0] += 2 * recipCoeff * force * kx ;
            forces[n][1] += 2 * recipCoeff * force * ky ;
            forces[n][2] += 2 * recipCoeff * force * kz ;
317
          }
318

319
          recipEnergy       = recipCoeff * ak * (cs * cs + ss * ss);
320
          totalRecipEnergy += recipEnergy;
321

322
          if (totalEnergy)
323
324
             *totalEnergy += recipEnergy;

325
326
          if (energyByAtom)
             for (int n = 0; n < numberOfAtoms; n++)
327
328
               energyByAtom[n] += recipEnergy;

329
330
331
332
333
          lowrz = 1 - numRz;
        }
        lowry = 1 - numRy;
      }
    }
334
  }
335
336
337
338
339

// **************************************************************************************
// SHORT-RANGE ENERGY AND FORCES
// **************************************************************************************

340
341
    if (!includeDirect)
        return;
peastman's avatar
peastman committed
342
343
    double totalVdwEnergy            = 0.0f;
    double totalRealSpaceEwaldEnergy = 0.0f;
344

345
346
347
348
    for (int i = 0; i < (int) neighborList->size(); i++) {
       OpenMM::AtomPair pair = (*neighborList)[i];
       int ii = pair.first;
       int jj = pair.second;
349

peastman's avatar
peastman committed
350
       double deltaR[2][ReferenceForce::LastDeltaRIndex];
351
       ReferenceForce::getDeltaRPeriodic(atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0]);
peastman's avatar
peastman committed
352
353
354
       double r         = deltaR[0][ReferenceForce::RIndex];
       double inverseR  = 1.0/(deltaR[0][ReferenceForce::RIndex]);
       double switchValue = 1, switchDeriv = 0;
355
       if (useSwitch && r > switchingDistance) {
peastman's avatar
peastman committed
356
           double t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
357
358
359
           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
360
       double alphaR    = alphaEwald * r;
361
362


peastman's avatar
peastman committed
363
364
       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);
365

peastman's avatar
peastman committed
366
367
368
369
370
371
372
       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;
373
374
375
376
       if (useSwitch) {
           dEdR -= vdwEnergy*switchDeriv*inverseR;
           vdwEnergy *= switchValue;
       }
377

378
       // accumulate forces
379

380
       for (int kk = 0; kk < 3; kk++) {
peastman's avatar
peastman committed
381
382
383
          double force = dEdR*deltaR[0][kk];
          forces[ii][kk] += force;
          forces[jj][kk] -= force;
384
385
       }

386
       // accumulate energies
387

peastman's avatar
peastman committed
388
       realSpaceEwaldEnergy        = ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR);
389

390
391
       totalVdwEnergy             += vdwEnergy;
       totalRealSpaceEwaldEnergy  += realSpaceEwaldEnergy;
392

393
        if (energyByAtom) {
394
395
396
           energyByAtom[ii] += realSpaceEwaldEnergy + vdwEnergy;
           energyByAtom[jj] += realSpaceEwaldEnergy + vdwEnergy;
        }
397

398
    }
399

400
    if (totalEnergy)
401
402
        *totalEnergy += totalRealSpaceEwaldEnergy + totalVdwEnergy;

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

peastman's avatar
peastman committed
405
    double totalExclusionEnergy = 0.0f;
406
    const double TWO_OVER_SQRT_PI = 2/sqrt(PI_M);
Peter Eastman's avatar
Peter Eastman committed
407
    for (int i = 0; i < numberOfAtoms; i++)
408
409
        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter) {
            if (*iter > i) {
Peter Eastman's avatar
Peter Eastman committed
410
               int ii = i;
411
               int jj = *iter;
Peter Eastman's avatar
Peter Eastman committed
412

peastman's avatar
peastman committed
413
               double deltaR[2][ReferenceForce::LastDeltaRIndex];
414
               ReferenceForce::getDeltaR(atomCoordinates[jj], atomCoordinates[ii], deltaR[0]);
peastman's avatar
peastman committed
415
416
417
               double r         = deltaR[0][ReferenceForce::RIndex];
               double inverseR  = 1.0/(deltaR[0][ReferenceForce::RIndex]);
               double alphaR    = alphaEwald * r;
418
               if (erf(alphaR) > 1e-6) {
peastman's avatar
peastman committed
419
420
                   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);
Peter Eastman's avatar
Peter Eastman committed
421

422
                   // accumulate forces
Peter Eastman's avatar
Peter Eastman committed
423

424
                   for (int kk = 0; kk < 3; kk++) {
peastman's avatar
peastman committed
425
                      double force  = dEdR*deltaR[0][kk];
426
427
428
                      forces[ii][kk]   -= force;
                      forces[jj][kk]   += force;
                   }
Peter Eastman's avatar
Peter Eastman committed
429

430
                   // accumulate energies
Peter Eastman's avatar
Peter Eastman committed
431

peastman's avatar
peastman committed
432
                   realSpaceEwaldEnergy = ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erf(alphaR);
433
434
               }
               else {
peastman's avatar
peastman committed
435
                   realSpaceEwaldEnergy = alphaEwald*TWO_OVER_SQRT_PI*ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex];
436
               }
Peter Eastman's avatar
Peter Eastman committed
437

438
439
440
441
               totalExclusionEnergy += realSpaceEwaldEnergy;
               if (energyByAtom) {
                   energyByAtom[ii] -= realSpaceEwaldEnergy;
                   energyByAtom[jj] -= realSpaceEwaldEnergy;
442
               }
Peter Eastman's avatar
Peter Eastman committed
443
            }
444
        }
Peter Eastman's avatar
Peter Eastman committed
445

446
    if (totalEnergy)
447
        *totalEnergy -= totalExclusionEnergy;
448
449
450
}


451
452
453
454
455
456
457
/**---------------------------------------------------------------------------------------

   Calculate LJ Coulomb pair ixn

   @param numberOfAtoms    number of atoms
   @param atomCoordinates  atom coordinates
   @param atomParameters   atom parameters                             atomParameters[atomIndex][paramterIndex]
458
459
   @param exclusions       atom exclusion indices
                           exclusions[atomIndex] contains the list of exclusions for that atom
460
461
462
463
   @param fixedParameters  non atom parameters (not currently used)
   @param forces           force array (forces added)
   @param energyByAtom     atom energy
   @param totalEnergy      total energy
464
465
   @param includeDirect      true if direct space interactions should be included
   @param includeReciprocal  true if reciprocal space interactions should be included
466
467

   --------------------------------------------------------------------------------------- */
468

peastman's avatar
peastman committed
469
470
471
472
void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<Vec3>& atomCoordinates,
                                             double** atomParameters, vector<set<int> >& exclusions,
                                             double* fixedParameters, vector<Vec3>& forces,
                                             double* energyByAtom, double* totalEnergy, bool includeDirect, bool includeReciprocal) const {
473

474
475
476
477
478
479
480
   if (ewald || pme) {
       calculateEwaldIxn(numberOfAtoms, atomCoordinates, atomParameters, exclusions, fixedParameters, forces, energyByAtom,
               totalEnergy, includeDirect, includeReciprocal);
       return;
   }
   if (!includeDirect)
       return;
481
   if (cutoff) {
482
       for (int i = 0; i < (int) neighborList->size(); i++) {
483
484
485
486
487
           OpenMM::AtomPair pair = (*neighborList)[i];
           calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
       }
   }
   else {
488
       for (int ii = 0; ii < numberOfAtoms; ii++) {
489
          // loop over atom pairs
490

491
          for (int jj = ii+1; jj < numberOfAtoms; jj++)
492
493
              if (exclusions[jj].find(ii) == exclusions[jj].end())
                  calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
494
495
496
       }
   }
}
497

498
  /**---------------------------------------------------------------------------------------
499

500
     Calculate LJ Coulomb pair ixn between two atoms
501

502
503
504
505
506
507
508
     @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 energyByAtom     atom energy
     @param totalEnergy      total energy
509

510
     --------------------------------------------------------------------------------------- */
511

peastman's avatar
peastman committed
512
513
514
void ReferenceLJCoulombIxn::calculateOneIxn(int ii, int jj, vector<Vec3>& atomCoordinates,
                        double** atomParameters, vector<Vec3>& forces,
                        double* energyByAtom, double* totalEnergy) const {
515

peastman's avatar
peastman committed
516
    double deltaR[2][ReferenceForce::LastDeltaRIndex];
517
518
519
520

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

    if (periodic)
521
        ReferenceForce::getDeltaRPeriodic(atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0]);
522
    else
523
        ReferenceForce::getDeltaR(atomCoordinates[jj], atomCoordinates[ii], deltaR[0]);
524

peastman's avatar
peastman committed
525
526
527
    double r2        = deltaR[0][ReferenceForce::R2Index];
    double inverseR  = 1.0/(deltaR[0][ReferenceForce::RIndex]);
    double switchValue = 1, switchDeriv = 0;
528
    if (useSwitch) {
peastman's avatar
peastman committed
529
        double r = deltaR[0][ReferenceForce::RIndex];
530
        if (r > switchingDistance) {
peastman's avatar
peastman committed
531
            double t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
532
533
534
535
            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
536
537
538
539
    double sig       = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
    double sig2      = inverseR*sig;
           sig2     *= sig2;
    double sig6      = sig2*sig2*sig2;
540

peastman's avatar
peastman committed
541
542
    double eps       = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
    double dEdR      = switchValue*eps*(12.0*sig6 - 6.0)*sig6;
543
    if (cutoff)
peastman's avatar
peastman committed
544
        dEdR += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
545
    else
peastman's avatar
peastman committed
546
        dEdR += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
547
    dEdR     *= inverseR*inverseR;
peastman's avatar
peastman committed
548
    double energy = eps*(sig6-1.0)*sig6;
549
550
551
552
553
    if (useSwitch) {
        dEdR -= energy*switchDeriv*inverseR;
        energy *= switchValue;
    }
    if (cutoff)
peastman's avatar
peastman committed
554
        energy += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf);
555
    else
peastman's avatar
peastman committed
556
        energy += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
557

558
559
    // accumulate forces

560
    for (int kk = 0; kk < 3; kk++) {
peastman's avatar
peastman committed
561
562
563
       double force = dEdR*deltaR[0][kk];
       forces[ii][kk] += force;
       forces[jj][kk] -= force;
564
565
566
567
    }

    // accumulate energies

568
    if (totalEnergy)
569
       *totalEnergy += energy;
570
    if (energyByAtom) {
571
572
       energyByAtom[ii] += energy;
       energyByAtom[jj] += energy;
573
574
    }
  }
575