CpuNonbondedForce.cpp 21.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
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
    if (pme) {
Peter Eastman's avatar
Peter Eastman committed
258
        ReferencePME pme(alphaEwald, numberOfAtoms, meshDim, 5, 1);
peastman's avatar
peastman committed
259
        vector<double> charges(numberOfAtoms);
peastman's avatar
peastman committed
260
261
        for (int i = 0; i < numberOfAtoms; i++)
            charges[i] = posq[4*i+3];
peastman's avatar
peastman committed
262
        double recipEnergy = 0.0;
Peter Eastman's avatar
Peter Eastman committed
263
        pme.exec(atomCoordinates, forces, charges, periodicBoxVectors, recipEnergy);
peastman's avatar
peastman committed
264
265
        if (totalEnergy)
            *totalEnergy += recipEnergy;
266
267
268

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

Peter Eastman's avatar
Peter Eastman committed
271
            vector<Vec3> dpmeforces;
272
            for (int i = 0; i < numberOfAtoms; i++){
peastman's avatar
peastman committed
273
274
                charges[i] = C6params[i];
                dpmeforces.push_back(Vec3());
275
            }
peastman's avatar
peastman committed
276
            double recipDispersionEnergy = 0.0;
Peter Eastman's avatar
Peter Eastman committed
277
            ljpme.exec_dpme(atomCoordinates, dpmeforces, charges, periodicBoxVectors, recipDispersionEnergy);
278
            for (int i = 0; i < numberOfAtoms; i++){
peastman's avatar
peastman committed
279
280
281
                forces[i][0] += dpmeforces[i][0];
                forces[i][1] += dpmeforces[i][1];
                forces[i][2] += dpmeforces[i][2];
282
283
284
285
286
            }
            if (totalEnergy)
                *totalEnergy += recipDispersionEnergy;
        }

peastman's avatar
peastman committed
287
    }
288
289
290

    // Ewald method

peastman's avatar
peastman committed
291
    else if (ewald) {
292

peastman's avatar
peastman committed
293
        // setup reciprocal box
294

peastman's avatar
peastman committed
295
        float recipBoxSize[3] = {(float) (TWO_PI/periodicBoxVectors[0][0]), (float) (TWO_PI/periodicBoxVectors[1][1]), (float) (TWO_PI/periodicBoxVectors[2][2])};
296
297


peastman's avatar
peastman committed
298
        // setup K-vectors
299

300
#define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
peastman's avatar
peastman committed
301
302
303
        vector<d_complex> eir(kmax*numberOfAtoms*3);
        vector<d_complex> tab_xy(numberOfAtoms);
        vector<d_complex> tab_qxyz(numberOfAtoms);
304

peastman's avatar
peastman committed
305
306
307
        for (int i = 0; (i < numberOfAtoms); i++) {
            float* pos = posq+4*i;
            for (int m = 0; (m < 3); m++)
308
                EIR(0, i, m) = d_complex(1,0);
309

peastman's avatar
peastman committed
310
            for (int m=0; (m<3); m++)
311
312
                EIR(1, i, m) = d_complex(cos(pos[m]*recipBoxSize[m]),
                                         sin(pos[m]*recipBoxSize[m]));
313

peastman's avatar
peastman committed
314
            for (int j=2; (j<kmax); j++)
315
316
                for (int m=0; (m<3); m++)
                    EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
317
318
        }

peastman's avatar
peastman committed
319
320
321
322
323
324
325
326
327
328
329
        // 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++)
330
                        tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
peastman's avatar
peastman committed
331
332
333
                }
                else {
                    for (int n = 0; n < numberOfAtoms; n++)
334
                        tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
peastman's avatar
peastman committed
335
336
337
338
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
                }
                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;
            }
371
372
        }
    }
peastman's avatar
peastman committed
373
}
374
375


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

408
409
410
411
412
413
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);
414
    float* forces = &(*threadForce)[threadIndex][0];
415
416
    fvec4 boxSize(periodicBoxVectors[0][0], periodicBoxVectors[1][1], periodicBoxVectors[2][2], 0);
    fvec4 invBoxSize(recipBoxSize[0], recipBoxSize[1], recipBoxSize[2], 0);
417
    if (ewald || pme || ljpme) {
418
        // Compute the interactions from the neighbor list.
419
        while (true) {
peastman's avatar
peastman committed
420
            int nextBlock = atomicCounter++;
421
422
423
424
            if (nextBlock >= neighborList->getNumBlocks())
                break;
            calculateBlockEwaldIxn(nextBlock, forces, energyPtr, boxSize, invBoxSize);
        }
425
426

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

peastman's avatar
peastman committed
428
429
        const int groupSize = max(1, numberOfAtoms/(10*numThreads));
        while (true) {
430
            int start = atomicCounter2.fetch_add(groupSize);
peastman's avatar
peastman committed
431
432
433
434
            if (start >= numberOfAtoms)
                break;
            int end = min(start+groupSize, numberOfAtoms);
            for (int i = start; i < end; i++) {
435
                fvec4 posI((float) atomCoordinates[i][0], (float) atomCoordinates[i][1], (float) atomCoordinates[i][2], 0.0f);
peastman's avatar
peastman committed
436
                float scaledChargeI = (float) (ONE_4PI_EPS0*posq[4*i+3]);
peastman's avatar
peastman committed
437
438
439
                for (int excluded : exclusions[i]) {
                    if (excluded > i) {
                        int j = excluded;
peastman's avatar
peastman committed
440
                        fvec4 deltaR;
peastman's avatar
peastman committed
441
                        fvec4 posJ((float) atomCoordinates[j][0], (float) atomCoordinates[j][1], (float) atomCoordinates[j][2], 0.0f);
peastman's avatar
peastman committed
442
                        float r2;
443
                        getDeltaR(posJ, posI, deltaR, r2, periodicExceptions, boxSize, invBoxSize);
peastman's avatar
peastman committed
444
445
446
447
448
449
450
                        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;
451
                            dEdR = dEdR * (erfAlphaR-TWO_OVER_SQRT_PI*alphaR*(float)exp(-alphaR*alphaR));
peastman's avatar
peastman committed
452
453
454
455
456
457
                            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;
                        }
458
                        else if (includeEnergy)
459
460
461
462
463
464
465
466
467
468
469
470
                            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);
                        }
471
                    }
472
                }
473
            }
474
        }
475
    }
476
    else {
477
        // Compute the interactions from the neighbor list.
478

479
        while (true) {
peastman's avatar
peastman committed
480
            int nextBlock = atomicCounter++;
481
482
483
484
            if (nextBlock >= neighborList->getNumBlocks())
                break;
            calculateBlockIxn(nextBlock, forces, energyPtr, boxSize, invBoxSize);
        }
485
    }
486
}
487

488
489
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;
490
    if (periodic) {
491
        if (triclinic) {
492
493
494
            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);
495
496
497
498
499
        }
        else {
            fvec4 base = round(deltaR*invBoxSize)*boxSize;
            deltaR = deltaR-base;
        }
500
    }
501
    r2 = dot3(deltaR, deltaR);
502
}
peastman's avatar
peastman committed
503

504
float CpuNonbondedForce::erfcApprox(float x) {
505
506
507
508
509
    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
510
}
peastman's avatar
peastman committed
511

512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
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];
}