CpuNonbondedForce.cpp 21.5 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
27
28
29
30
 * 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 <complex>

#include "SimTKOpenMMUtilities.h"
#include "CpuNonbondedForce.h"
#include "ReferenceForce.h"
#include "ReferencePME.h"
31
#include <algorithm>
32
#include <iostream>
33
34
35
36
37
38

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

using namespace std;
peastman's avatar
peastman committed
39
using namespace OpenMM;
40

peastman's avatar
peastman committed
41
const float CpuNonbondedForce::TWO_OVER_SQRT_PI = (float) (2/sqrt(PI_M));
42
const int CpuNonbondedForce::NUM_TABLE_POINTS = 2048;
43

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

   CpuNonbondedForce constructor

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

50
51
52
CpuNonbondedForce::CpuNonbondedForce(const CpuNeighborList& neighbors) : neighborList(&neighbors), cutoff(false), useSwitch(false), periodic(false),
        periodicExceptions(false), ewald(false), pme(false), ljpme(false), tableIsValid(false), expTableIsValid(false), cutoffDistance(0.0f),
        alphaDispersionEwald(0.0f), alphaEwald(0.0f) {
53
54
}

55
56
57
CpuNonbondedForce::~CpuNonbondedForce() {
}

58
/**---------------------------------------------------------------------------------------
59

60
   Set the force to use a cutoff.
61

62
63
   @param distance            the cutoff distance
   @param solventDielectric   the dielectric constant of the bulk solvent
64
65
66

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

67
void CpuNonbondedForce::setUseCutoff(float distance, float solventDielectric) {
68
69
    if (distance != cutoffDistance)
        tableIsValid = false;
70
71
    cutoff = true;
    cutoffDistance = distance;
72
    inverseRcut6 = pow(cutoffDistance, -6);
73
74
    krf = pow(cutoffDistance, -3.0f)*(solventDielectric-1.0)/(2.0*solventDielectric+1.0);
    crf = (1.0/cutoffDistance)*(3.0*solventDielectric)/(2.0*solventDielectric+1.0);
75
76
77
78
79
80
81
82
83
84
85
    if(alphaDispersionEwald != 0.0f){
        // We set this here, in case setUseCutoff is called after the dispersion alpha is set.
        double dalphaR = alphaDispersionEwald*cutoffDistance;
        double dar2 = dalphaR * dalphaR;
        double dar4 = dar2*dar2;
        double dar6 = dar4*dar2;
        double expterm = EXP(-dar2);
        inverseRcut6Expterm  = inverseRcut6*(1.0 - expterm * (1.0 + dar2 + 0.5*dar4));
    }

}
86
87
88
89
90
91
92
93
94
95
96
97
98
99

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

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

   @param distance            the switching distance

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

void CpuNonbondedForce::setUseSwitchingFunction(float distance) {
    useSwitch = true;
    switchingDistance = distance;
}

100
/**---------------------------------------------------------------------------------------
101
102
103
104
105

     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.

106
     @param periodicBoxVectors    the vectors defining the periodic box
107
108
109

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

peastman's avatar
peastman committed
110
void CpuNonbondedForce::setPeriodic(Vec3* periodicBoxVectors) {
111
112

    assert(cutoff);
113
114
115
    assert(periodicBoxVectors[0][0] >= 2.0*cutoffDistance);
    assert(periodicBoxVectors[1][1] >= 2.0*cutoffDistance);
    assert(periodicBoxVectors[2][2] >= 2.0*cutoffDistance);
116
    periodic = true;
117
118
119
120
121
122
123
124
125
126
127
    this->periodicBoxVectors[0] = periodicBoxVectors[0];
    this->periodicBoxVectors[1] = periodicBoxVectors[1];
    this->periodicBoxVectors[2] = periodicBoxVectors[2];
    recipBoxSize[0] = (float) (1.0/periodicBoxVectors[0][0]);
    recipBoxSize[1] = (float) (1.0/periodicBoxVectors[1][1]);
    recipBoxSize[2] = (float) (1.0/periodicBoxVectors[2][2]);
    periodicBoxVec4.resize(3);
    periodicBoxVec4[0] = fvec4(periodicBoxVectors[0][0], periodicBoxVectors[0][1], periodicBoxVectors[0][2], 0);
    periodicBoxVec4[1] = fvec4(periodicBoxVectors[1][0], periodicBoxVectors[1][1], periodicBoxVectors[1][2], 0);
    periodicBoxVec4[2] = fvec4(periodicBoxVectors[2][0], periodicBoxVectors[2][1], periodicBoxVectors[2][2], 0);
    triclinic = (periodicBoxVectors[0][1] != 0.0 || periodicBoxVectors[0][2] != 0.0 ||
128
129
130
            periodicBoxVectors[1][0] != 0.0 || periodicBoxVectors[1][2] != 0.0 ||
            periodicBoxVectors[2][0] != 0.0 || periodicBoxVectors[2][1] != 0.0);
}
131

132
/**---------------------------------------------------------------------------------------
133
134
135
136
137
138
139
140
141
142

     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

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

143
144
145
146
147
148
149
150
151
152
void CpuNonbondedForce::setUseEwald(float alpha, int kmaxx, int kmaxy, int kmaxz) {
    if (alpha != alphaEwald)
        tableIsValid = false;
    alphaEwald = alpha;
    numRx = kmaxx;
    numRy = kmaxy;
    numRz = kmaxz;
    ewald = true;
    tabulateEwaldScaleFactor();
}
153

154
/**---------------------------------------------------------------------------------------
155
156
157
158
159
160
161
162

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

     @param alpha  the Ewald separation parameter
     @param gridSize the dimensions of the mesh

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

163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
void CpuNonbondedForce::setUsePME(float alpha, int meshSize[3]) {
    if (alpha != alphaEwald)
        tableIsValid = false;
    alphaEwald = alpha;
    meshDim[0] = meshSize[0];
    meshDim[1] = meshSize[1];
    meshDim[2] = meshSize[2];
    pme = true;
    tabulateEwaldScaleFactor();
}


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

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

     @param alpha  the Ewald separation parameter
     @param gridSize the dimensions of the mesh

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

void CpuNonbondedForce::setUseLJPME(float alpha, int meshSize[3]) {
    if (alpha != alphaDispersionEwald)
        expTableIsValid = false;
    alphaDispersionEwald = alpha;
    dispersionMeshDim[0] = meshSize[0];
    dispersionMeshDim[1] = meshSize[1];
    dispersionMeshDim[2] = meshSize[2];
    ljpme = true;
    tabulateExpTerms();
    if(cutoffDistance != 0.0f){
        // We set this here, in case setUseLJPME is called after the cutoff is set
        double dalphaR = alphaDispersionEwald*cutoffDistance;
        double dar2 = dalphaR * dalphaR;
        double dar4 = dar2*dar2;
        double dar6 = dar4*dar2;
        double expterm = EXP(-dar2);
        inverseRcut6Expterm  = inverseRcut6*(1.0 - expterm * (1.0 + dar2 + 0.5*dar4));
    }
}

204
205
206
void CpuNonbondedForce::setPeriodicExceptions(bool periodic) {
    periodicExceptions = periodic;
}
207
208

void CpuNonbondedForce::tabulateEwaldScaleFactor() {
209
210
211
    if (tableIsValid)
        return;
    tableIsValid = true;
peastman's avatar
peastman committed
212
    ewaldDX = cutoffDistance/NUM_TABLE_POINTS;
peastman's avatar
peastman committed
213
    ewaldDXInv = 1.0f/ewaldDX;
214
215
    erfcDXInv = 1.0f/(ewaldDX*alphaEwald);
    erfcTable.resize(NUM_TABLE_POINTS+4);
peastman's avatar
peastman committed
216
217
    ewaldScaleTable.resize(NUM_TABLE_POINTS+4);
    for (int i = 0; i < NUM_TABLE_POINTS+4; i++) {
218
        double r = i*ewaldDX;
peastman's avatar
peastman committed
219
        double alphaR = alphaEwald*r;
220
221
        erfcTable[i] = erfc(alphaR);
        ewaldScaleTable[i] = erfcTable[i] + TWO_OVER_SQRT_PI*alphaR*exp(-alphaR*alphaR);
peastman's avatar
peastman committed
222
223
    }
}
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244

void CpuNonbondedForce::tabulateExpTerms() {
    if (expTableIsValid)
        return;
    expTableIsValid = true;
    exptermsDX = cutoffDistance/NUM_TABLE_POINTS;
    exptermsDXInv = 1.0f/exptermsDX;
    exptermsTable.resize(NUM_TABLE_POINTS+4);
    dExptermsTable.resize(NUM_TABLE_POINTS+4);
    for (int i = 0; i < NUM_TABLE_POINTS+4; i++) {
        double r = i*ewaldDX;
        double dalphaR = alphaDispersionEwald*r;
        double dar2 = dalphaR * dalphaR;
        double dar4 = dar2*dar2;
        double dar6 = dar4*dar2;
        double expterm = EXP(-dar2);
        exptermsTable[i]  = (1.0 - expterm * (1.0 + dar2 + 0.5*dar4));
        dExptermsTable[i] = (1.0 - expterm * (1.0 + dar2 + 0.5*dar4 + dar6/6.0));
    }
}

peastman's avatar
peastman committed
245
void CpuNonbondedForce::calculateReciprocalIxn(int numberOfAtoms, float* posq, const vector<Vec3>& atomCoordinates,
246
                                               const vector<pair<float, float> >& atomParameters, const vector<float> &C6params, const vector<set<int> >& exclusions,
peastman's avatar
peastman committed
247
                                               vector<Vec3>& forces, double* totalEnergy) const {
248
249
250
251
    typedef std::complex<float> d_complex;

    static const float epsilon     =  1.0;

252
    int kmax                       = (ewald ? max(numRx, max(numRy,numRz)) : 0);
peastman's avatar
peastman committed
253
    float factorEwald              = -1 / (4*alphaEwald*alphaEwald);
254
    float TWO_PI                   = 2.0 * PI_M;
255
    float recipCoeff               = (float)(ONE_4PI_EPS0*4*PI_M/(periodicBoxVectors[0][0] * periodicBoxVectors[1][1] * periodicBoxVectors[2][2]) /epsilon);
256

peastman's avatar
peastman committed
257
258
259
    if (pme) {
        pme_t pmedata;
        pme_init(&pmedata, alphaEwald, numberOfAtoms, meshDim, 5, 1);
peastman's avatar
peastman committed
260
        vector<double> charges(numberOfAtoms);
peastman's avatar
peastman committed
261
262
        for (int i = 0; i < numberOfAtoms; i++)
            charges[i] = posq[4*i+3];
peastman's avatar
peastman committed
263
        double recipEnergy = 0.0;
264
        pme_exec(pmedata, atomCoordinates, forces, charges, periodicBoxVectors, &recipEnergy);
peastman's avatar
peastman committed
265
266
        if (totalEnergy)
            *totalEnergy += recipEnergy;
267
        pme_destroy(pmedata);
268
269
270
271
272

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

peastman's avatar
peastman committed
273
            std::vector<Vec3> dpmeforces;
274
            for (int i = 0; i < numberOfAtoms; i++){
peastman's avatar
peastman committed
275
276
                charges[i] = C6params[i];
                dpmeforces.push_back(Vec3());
277
            }
peastman's avatar
peastman committed
278
            double recipDispersionEnergy = 0.0;
279
280
            pme_exec_dpme(pmedata,atomCoordinates,dpmeforces,charges,periodicBoxVectors,&recipDispersionEnergy);
            for (int i = 0; i < numberOfAtoms; i++){
peastman's avatar
peastman committed
281
282
283
                forces[i][0] += dpmeforces[i][0];
                forces[i][1] += dpmeforces[i][1];
                forces[i][2] += dpmeforces[i][2];
284
285
286
287
288
289
290
            }
            if (totalEnergy)
                *totalEnergy += recipDispersionEnergy;

            pme_destroy(pmedata);
        }

peastman's avatar
peastman committed
291
    }
292
293
294

    // Ewald method

peastman's avatar
peastman committed
295
    else if (ewald) {
296

peastman's avatar
peastman committed
297
        // setup reciprocal box
298

peastman's avatar
peastman committed
299
        float recipBoxSize[3] = {(float) (TWO_PI/periodicBoxVectors[0][0]), (float) (TWO_PI/periodicBoxVectors[1][1]), (float) (TWO_PI/periodicBoxVectors[2][2])};
300
301


peastman's avatar
peastman committed
302
        // setup K-vectors
303

304
#define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
peastman's avatar
peastman committed
305
306
307
        vector<d_complex> eir(kmax*numberOfAtoms*3);
        vector<d_complex> tab_xy(numberOfAtoms);
        vector<d_complex> tab_qxyz(numberOfAtoms);
308

peastman's avatar
peastman committed
309
310
311
        for (int i = 0; (i < numberOfAtoms); i++) {
            float* pos = posq+4*i;
            for (int m = 0; (m < 3); m++)
312
                EIR(0, i, m) = d_complex(1,0);
313

peastman's avatar
peastman committed
314
            for (int m=0; (m<3); m++)
315
316
                EIR(1, i, m) = d_complex(cos(pos[m]*recipBoxSize[m]),
                                         sin(pos[m]*recipBoxSize[m]));
317

peastman's avatar
peastman committed
318
            for (int j=2; (j<kmax); j++)
319
320
                for (int m=0; (m<3); m++)
                    EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
321
322
        }

peastman's avatar
peastman committed
323
324
325
326
327
328
329
330
331
332
333
        // calculate reciprocal space energy and forces

        int lowry = 0;
        int lowrz = 1;

        for (int rx = 0; rx < numRx; rx++) {
            float kx = rx * recipBoxSize[0];
            for (int ry = lowry; ry < numRy; ry++) {
                float ky = ry * recipBoxSize[1];
                if (ry >= 0) {
                    for (int n = 0; n < numberOfAtoms; n++)
334
                        tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
peastman's avatar
peastman committed
335
336
337
                }
                else {
                    for (int n = 0; n < numberOfAtoms; n++)
338
                        tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
peastman's avatar
peastman committed
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
                }
                for (int rz = lowrz; rz < numRz; rz++) {
                    if (rz >= 0) {
                        for (int n = 0; n < numberOfAtoms; n++)
                            tab_qxyz[n] = posq[4*n+3] * (tab_xy[n] * EIR(rz, n, 2));
                    }
                    else {
                        for (int n = 0; n < numberOfAtoms; n++)
                            tab_qxyz[n] = posq[4*n+3] * (tab_xy[n] * conj(EIR(-rz, n, 2)));
                    }
                    float cs = 0.0f;
                    float ss = 0.0f;

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

                    float kz = rz * recipBoxSize[2];
                    float k2 = kx * kx + ky * ky + kz * kz;
                    float ak = exp(k2*factorEwald) / k2;

                    for (int n = 0; n < numberOfAtoms; n++) {
                        float force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real());
                        forces[n][0] += 2 * recipCoeff * force * kx;
                        forces[n][1] += 2 * recipCoeff * force * ky;
                        forces[n][2] += 2 * recipCoeff * force * kz;
                    }

                    if (totalEnergy)
                        *totalEnergy += recipCoeff * ak * (cs * cs + ss * ss);

                    lowrz = 1 - numRz;
                }
                lowry = 1 - numRy;
            }
375
376
        }
    }
peastman's avatar
peastman committed
377
}
378
379


peastman's avatar
peastman committed
380
void CpuNonbondedForce::calculateDirectIxn(int numberOfAtoms, float* posq, const vector<Vec3>& atomCoordinates, const vector<pair<float, float> >& atomParameters,
381
                                           const vector<float>& C6params, const vector<set<int> >& exclusions, vector<AlignedArray<float> >& threadForce, double* totalEnergy, ThreadPool& threads) {
382
383
    // Record the parameters for the threads.
    
peastman's avatar
peastman committed
384
    this->numberOfAtoms = numberOfAtoms;
385
    this->posq = posq;
386
    this->atomCoordinates = &atomCoordinates[0];
peastman's avatar
peastman committed
387
    this->atomParameters = &atomParameters[0];
388
    this->C6params = &C6params[0];
peastman's avatar
peastman committed
389
    this->exclusions = &exclusions[0];
390
    this->threadForce = &threadForce;
391
    includeEnergy = (totalEnergy != NULL);
392
    threadEnergy.resize(threads.getNumThreads());
peastman's avatar
peastman committed
393
    atomicCounter = 0;
394
395
396
    
    // Signal the threads to start running and wait for them to finish.
    
peastman's avatar
peastman committed
397
    threads.execute([&] (ThreadPool& threads, int threadIndex) { threadComputeDirect(threads, threadIndex); });
398
    threads.waitForThreads();
399
    
peastman's avatar
peastman committed
400
401
402
    // Signal the threads to subtract the exclusions.
    
    if (ewald || pme) {
peastman's avatar
peastman committed
403
        atomicCounter = 0;
peastman's avatar
peastman committed
404
405
406
407
        threads.resumeThreads();
        threads.waitForThreads();
    }
    
408
    // Combine the energies from all the threads.
409
    
410
411
412
413
414
    if (totalEnergy != NULL) {
        double directEnergy = 0;
        int numThreads = threads.getNumThreads();
        for (int i = 0; i < numThreads; i++)
            directEnergy += threadEnergy[i];
415
        *totalEnergy += directEnergy;
416
    }
417
418
}

419
420
421
422
423
424
void CpuNonbondedForce::threadComputeDirect(ThreadPool& threads, int threadIndex) {
    // Compute this thread's subset of interactions.

    int numThreads = threads.getNumThreads();
    threadEnergy[threadIndex] = 0;
    double* energyPtr = (includeEnergy ? &threadEnergy[threadIndex] : NULL);
425
    float* forces = &(*threadForce)[threadIndex][0];
426
427
    fvec4 boxSize(periodicBoxVectors[0][0], periodicBoxVectors[1][1], periodicBoxVectors[2][2], 0);
    fvec4 invBoxSize(recipBoxSize[0], recipBoxSize[1], recipBoxSize[2], 0);
428
    if (ewald || pme || ljpme) {
429
        // Compute the interactions from the neighbor list.
430
        while (true) {
peastman's avatar
peastman committed
431
            int nextBlock = atomicCounter++;
432
433
434
435
            if (nextBlock >= neighborList->getNumBlocks())
                break;
            calculateBlockEwaldIxn(nextBlock, forces, energyPtr, boxSize, invBoxSize);
        }
436
437

        // Now subtract off the exclusions, since they were implicitly included in the reciprocal space sum.
438

peastman's avatar
peastman committed
439
440
441
        threads.syncThreads();
        const int groupSize = max(1, numberOfAtoms/(10*numThreads));
        while (true) {
peastman's avatar
peastman committed
442
            int start = atomicCounter.fetch_add(groupSize);
peastman's avatar
peastman committed
443
444
445
446
            if (start >= numberOfAtoms)
                break;
            int end = min(start+groupSize, numberOfAtoms);
            for (int i = start; i < end; i++) {
447
                fvec4 posI((float) atomCoordinates[i][0], (float) atomCoordinates[i][1], (float) atomCoordinates[i][2], 0.0f);
peastman's avatar
peastman committed
448
                float scaledChargeI = (float) (ONE_4PI_EPS0*posq[4*i+3]);
peastman's avatar
peastman committed
449
450
451
                for (int excluded : exclusions[i]) {
                    if (excluded > i) {
                        int j = excluded;
peastman's avatar
peastman committed
452
                        fvec4 deltaR;
peastman's avatar
peastman committed
453
                        fvec4 posJ((float) atomCoordinates[j][0], (float) atomCoordinates[j][1], (float) atomCoordinates[j][2], 0.0f);
peastman's avatar
peastman committed
454
                        float r2;
455
                        getDeltaR(posJ, posI, deltaR, r2, periodicExceptions, boxSize, invBoxSize);
peastman's avatar
peastman committed
456
457
458
459
460
461
462
                        float r = sqrtf(r2);
                        float alphaR = alphaEwald*r;
                        float erfAlphaR = erf(alphaR);
                        if (erfAlphaR > 1e-6f) {
                            float inverseR = 1/r;
                            float chargeProdOverR = scaledChargeI*posq[4*j+3]*inverseR;
                            float dEdR = chargeProdOverR*inverseR*inverseR;
463
                            dEdR = dEdR * (erfAlphaR-TWO_OVER_SQRT_PI*alphaR*(float)exp(-alphaR*alphaR));
peastman's avatar
peastman committed
464
465
466
467
468
469
                            fvec4 result = deltaR*dEdR;
                            (fvec4(forces+4*i)-result).store(forces+4*i);
                            (fvec4(forces+4*j)+result).store(forces+4*j);
                            if (includeEnergy)
                                threadEnergy[threadIndex] -= chargeProdOverR*erfAlphaR;
                        }
470
                        else if (includeEnergy)
471
472
473
474
475
476
477
478
479
480
481
482
                            threadEnergy[threadIndex] -= alphaEwald*TWO_OVER_SQRT_PI*scaledChargeI*posq[4*j+3];
                        if (ljpme) {
                            float C6ij = C6params[i]*C6params[j];
                            float inverseR2 = 1.0f/r2;
                            float emult = C6ij*inverseR2*inverseR2*inverseR2*exptermsApprox(r);
                            if(includeEnergy)
                                threadEnergy[threadIndex] += emult;
                            float dEdR = -6.0f*C6ij*inverseR2*inverseR2*inverseR2*inverseR2*dExptermsApprox(r);
                            fvec4 result = deltaR*dEdR;
                            (fvec4(forces+4*i)-result).store(forces+4*i);
                            (fvec4(forces+4*j)+result).store(forces+4*j);
                        }
483
                    }
484
                }
485
            }
486
        }
487
    }
488
    else {
489
        // Compute the interactions from the neighbor list.
490

491
        while (true) {
peastman's avatar
peastman committed
492
            int nextBlock = atomicCounter++;
493
494
495
496
            if (nextBlock >= neighborList->getNumBlocks())
                break;
            calculateBlockIxn(nextBlock, forces, energyPtr, boxSize, invBoxSize);
        }
497
    }
498
}
499

500
501
void CpuNonbondedForce::getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, bool periodic, const fvec4& boxSize, const fvec4& invBoxSize) const {
    deltaR = posJ-posI;
502
    if (periodic) {
503
        if (triclinic) {
504
505
506
            deltaR -= periodicBoxVec4[2]*floorf(deltaR[2]*recipBoxSize[2]+0.5f);
            deltaR -= periodicBoxVec4[1]*floorf(deltaR[1]*recipBoxSize[1]+0.5f);
            deltaR -= periodicBoxVec4[0]*floorf(deltaR[0]*recipBoxSize[0]+0.5f);
507
508
509
510
511
        }
        else {
            fvec4 base = round(deltaR*invBoxSize)*boxSize;
            deltaR = deltaR-base;
        }
512
    }
513
    r2 = dot3(deltaR, deltaR);
514
}
peastman's avatar
peastman committed
515

516
float CpuNonbondedForce::erfcApprox(float x) {
517
518
519
520
521
    float x1 = x*erfcDXInv;
    int index = min((int) floor(x1), NUM_TABLE_POINTS);
    float coeff2 = x1-index;
    float coeff1 = 1.0f-coeff2;
    return coeff1*erfcTable[index] + coeff2*erfcTable[index+1];
peastman's avatar
peastman committed
522
}
peastman's avatar
peastman committed
523

524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
float CpuNonbondedForce::exptermsApprox(float x) {
    float x1 = x*exptermsDXInv;
    int index = min((int) floor(x1), NUM_TABLE_POINTS);
    float coeff2 = x1-index;
    float coeff1 = 1.0f-coeff2;
    return coeff1*exptermsTable[index] + coeff2*exptermsTable[index+1];
}

float CpuNonbondedForce::dExptermsApprox(float x) {
    float x1 = x*exptermsDXInv;
    int index = min((int) floor(x1), NUM_TABLE_POINTS);
    float coeff2 = x1-index;
    float coeff1 = 1.0f-coeff2;
    return coeff1*dExptermsTable[index] + coeff2*dExptermsTable[index+1];
}