CpuPmeKernels.cpp 38.7 KB
Newer Older
peastman's avatar
peastman committed
1
2
3
4
5
6
7
8
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
9
 * Portions copyright (c) 2013-2025 Stanford University and the Authors.      *
peastman's avatar
peastman committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
 * Authors: Peter Eastman                                                     *
 * Contributors:                                                              *
 *                                                                            *
 * 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.                                     *
 * -------------------------------------------------------------------------- */

32
33
34
#ifdef WIN32
  #define _USE_MATH_DEFINES // Needed to get M_PI
#endif
Peter Eastman's avatar
Peter Eastman committed
35
36
37
38
#ifdef _MSC_VER
  #define POCKETFFT_NO_VECTORS
#endif
#define POCKETFFT_CACHE_SIZE 4
39
#include "CpuPmeKernels.h"
40
#include "SimTKOpenMMRealType.h"
41
#include "ReferenceForce.h"
42
#include "openmm/internal/hardware.h"
43
#include "openmm/internal/vectorize.h"
44
#include "openmm/OpenMMException.h"
Peter Eastman's avatar
Peter Eastman committed
45
#include "pocketfft_hdronly.h"
46
#include <cmath>
47
#include <algorithm>
peastman's avatar
peastman committed
48
#include <cstring>
49
#include <sstream>
Robert McGibbon's avatar
Robert McGibbon committed
50
#include <cstdlib>
peastman's avatar
peastman committed
51
52
53
54
55
56

using namespace OpenMM;
using namespace std;

static const int PME_ORDER = 5;

57
58
bool CpuCalcDispersionPmeReciprocalForceKernel::hasInitializedThreads = false;
int CpuCalcDispersionPmeReciprocalForceKernel::numThreads = 0;
59

Peter Eastman's avatar
Peter Eastman committed
60
static void spreadCharge(float* posq, vector<float>& grid, int gridx, int gridy, int gridz, int numParticles, Vec3* periodicBoxVectors, Vec3* recipBoxVectors,
peastman's avatar
peastman committed
61
        atomic<int>& atomicCounter, const float epsilonFactor, int threadIndex, int numThreads, bool deterministic) {
62
    float temp[4];
peastman's avatar
peastman committed
63
64
65
66
67
    fvec4 boxSize((float) periodicBoxVectors[0][0], (float) periodicBoxVectors[1][1], (float) periodicBoxVectors[2][2], 0);
    fvec4 invBoxSize((float) recipBoxVectors[0][0], (float) recipBoxVectors[1][1], (float) recipBoxVectors[2][2], 0);
    fvec4 recipBoxVec0((float) recipBoxVectors[0][0], (float) recipBoxVectors[0][1], (float) recipBoxVectors[0][2], 0);
    fvec4 recipBoxVec1((float) recipBoxVectors[1][0], (float) recipBoxVectors[1][1], (float) recipBoxVectors[1][2], 0);
    fvec4 recipBoxVec2((float) recipBoxVectors[2][0], (float) recipBoxVectors[2][1], (float) recipBoxVectors[2][2], 0);
68
69
70
71
    fvec4 gridSize(gridx, gridy, gridz, 0);
    ivec4 gridSizeInt(gridx, gridy, gridz, 0);
    fvec4 one(1);
    fvec4 scale(1.0f/(PME_ORDER-1));
Robert McGibbon's avatar
Robert McGibbon committed
72
    float posInBox[4] = {0,0,0,0};
Peter Eastman's avatar
Peter Eastman committed
73
    memset(grid.data(), 0, sizeof(float)*gridx*gridy*gridz);
Robert McGibbon's avatar
Robert McGibbon committed
74

75
76
    const int groupSize = max(1, numParticles / (10 * numThreads));
    int start = groupSize * threadIndex;
77
    while (true) {
78
        if (!deterministic)
79
80
81
            start = atomicCounter.fetch_add(groupSize);

        if (start >= numParticles)
82
83
            break;

84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
        int end = min(start + groupSize, numParticles);
        for (int i = start; i < end; ++i) {
            // Find the position relative to the nearest grid point.

            fvec4 pos(&posq[4*i]);
            (pos-boxSize*floor(pos*invBoxSize)).store(posInBox);
            fvec4 t = posInBox[0]*recipBoxVec0 + posInBox[1]*recipBoxVec1 + posInBox[2]*recipBoxVec2;
            t = (t-floor(t))*gridSize;
            ivec4 ti = t;
            fvec4 dr = t-ti;
            ivec4 gridIndex = ti-(gridSizeInt&ti==gridSizeInt);

            // Compute the B-spline coefficients.

            fvec4 data[PME_ORDER];
            data[PME_ORDER-1] = 0.0f;
            data[1] = dr;
            data[0] = one-dr;
            for (int j = 3; j < PME_ORDER; j++) {
                fvec4 div(1.0f/(j-1));
                data[j-1] = div*dr*data[j-2];
                for (int k = 1; k < j-1; k++)
                    data[j-k-1] = div*((dr+k)*data[j-k-2]+(fvec4(j-k)-dr)*data[j-k-1]);
                data[0] = div*(one-dr)*data[0];
            }
            data[PME_ORDER-1] = scale*dr*data[PME_ORDER-2];
            for (int j = 1; j < (PME_ORDER-1); j++)
                data[PME_ORDER-j-1] = scale*((dr+j)*data[PME_ORDER-j-2]+(fvec4(PME_ORDER-j)-dr)*data[PME_ORDER-j-1]);
            data[0] = scale*(one-dr)*data[0];

            // Spread the charges.

            int gridIndexX = gridIndex[0];
            int gridIndexY = gridIndex[1];
            int gridIndexZ = gridIndex[2];
            if (gridIndexX < 0)
                return; // This happens when a simulation blows up and coordinates become NaN.
            int zindex[PME_ORDER];
            for (int j = 0; j < PME_ORDER; j++) {
                zindex[j] = gridIndexZ+j;
                zindex[j] -= (zindex[j] >= gridz ? gridz : 0);
            }
            float charge = epsilonFactor*posq[4*i+3];
            fvec4 zdata0to3(data[0][2], data[1][2], data[2][2], data[3][2]);
            float zdata4 = data[4][2];
            if (gridIndexZ+4 < gridz) {
                for (int ix = 0; ix < PME_ORDER; ix++) {
                    int xbase = gridIndexX+ix;
                    xbase -= (xbase >= gridx ? gridx : 0);
                    xbase = xbase*gridy*gridz;
                    float xdata = charge*data[ix][0];
                    for (int iy = 0; iy < PME_ORDER; iy++) {
                        int ybase = gridIndexY+iy;
                        ybase -= (ybase >= gridy ? gridy : 0);
                        ybase = xbase + ybase*gridz;
                        float multiplier = xdata*data[iy][1];
                        fvec4 add0to3 = zdata0to3*multiplier;
                        (fvec4(&grid[ybase+gridIndexZ])+add0to3).store(&grid[ybase+gridIndexZ]);
                        grid[ybase+zindex[4]] += multiplier*zdata4;
                    }
144
145
                }
            }
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
            else {
                for (int ix = 0; ix < PME_ORDER; ix++) {
                    int xbase = gridIndexX+ix;
                    xbase -= (xbase >= gridx ? gridx : 0);
                    xbase = xbase*gridy*gridz;
                    float xdata = charge*data[ix][0];
                    for (int iy = 0; iy < PME_ORDER; iy++) {
                        int ybase = gridIndexY+iy;
                        ybase -= (ybase >= gridy ? gridy : 0);
                        ybase = xbase + ybase*gridz;
                        float multiplier = xdata*data[iy][1];
                        fvec4 add0to3 = zdata0to3*multiplier;
                        add0to3.store(temp);
                        grid[ybase+zindex[0]] += temp[0];
                        grid[ybase+zindex[1]] += temp[1];
                        grid[ybase+zindex[2]] += temp[2];
                        grid[ybase+zindex[3]] += temp[3];
                        grid[ybase+zindex[4]] += multiplier*zdata4;
                    }
peastman's avatar
peastman committed
165
166
167
                }
            }
        }
168

169
        if (deterministic)
170
            start += groupSize * numThreads;
peastman's avatar
peastman committed
171
172
173
    }
}

174
#define FAST_ERFC 1
175
176
177
static void computeReciprocalDispersionEterm(int start, int end, int gridx, int gridy, int gridz, vector<float>& recipEterm, double alpha, vector<float>* bsplineModuli, Vec3* periodicBoxVectors, Vec3* recipBoxVectors) {
    const unsigned int zsize = gridz/2+1;
    const unsigned int yzsize = gridy*zsize;
178
    const float scaleFactor = (float)  -2.0f*M_PI*sqrtf(M_PI) / (6.0*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2]);
179
180
181
182
183

    float bfac = M_PI / alpha;
    float fac1 = 2.0f*M_PI*M_PI*M_PI*sqrtf(M_PI);
    float fac2 = alpha*alpha*alpha;
    float fac3 = -2.0f*alpha*M_PI*M_PI;
184
    float b, m, m3, expterm, erfcterm, t;
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205

    for (int kx = start; kx < end; kx++) {
        int mx = (kx < (gridx+1)/2) ? kx : kx-gridx;
        float mhx = mx*(float)recipBoxVectors[0][0];
        float bx = bsplineModuli[0][kx];
        for (int ky = 0; ky < gridy; ky++) {
            int my = (ky < (gridy+1)/2) ? ky : ky-gridy;
            float mhy = mx*(float)recipBoxVectors[1][0] + my*(float)recipBoxVectors[1][1];
            float mhx2y2 = mhx*mhx + mhy*mhy;
            float bxby = bx*bsplineModuli[1][ky];
            for (int kz = 0; kz < zsize; kz++) {
                int index = kx*yzsize + ky*zsize + kz;
                int mz = (kz < (gridz+1)/2) ? kz : kz-gridz;
                float mhz = mx*(float)recipBoxVectors[2][0] + my*(float)recipBoxVectors[2][1] + mz*(float)recipBoxVectors[2][2];
                float bz = bsplineModuli[2][kz];
                float m2 = mhx2y2 + mhz*mhz;
                float denom = scaleFactor/(bxby*bz);

                m = sqrtf(m2);
                m3 = m*m2;
                b = bfac*m;
206
207
208
209
210
211
212
213
214
                expterm = exp(-b*b);

#if FAST_ERFC
                // This approximation for erfc is from Abramowitz and Stegun (1964) p. 299.  They cite the following as
                // the original source: C. Hastings, Jr., Approximations for Digital Computers (1955).  It has a maximum
                // error of 1.5e-7.  Stolen by ACS from the CUDA platform's AMOEBA plugin.
                t = 1.0f/(1.0f+0.3275911f*b);
                erfcterm = (0.254829592f+(-0.284496736f+(1.421413741f+(-1.453152027f+1.061405429f*t)*t)*t)*t)*t*expterm;
#else
215
                erfcterm = erfc(b);
216
217
#endif
                recipEterm[index] = (fac1*erfcterm*m3 + expterm*(fac2 + fac3*m2)) * denom;
218
219
220
221
222
            }
        }
    }
}

peastman's avatar
peastman committed
223
static void computeReciprocalEterm(int start, int end, int gridx, int gridy, int gridz, vector<float>& recipEterm, double alpha, vector<float>* bsplineModuli, Vec3* periodicBoxVectors, Vec3* recipBoxVectors) {
224
225
    const unsigned int zsize = gridz/2+1;
    const unsigned int yzsize = gridy*zsize;
peastman's avatar
peastman committed
226
    const float scaleFactor = (float) (M_PI*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2]);
227
228
    const float recipExpFactor = (float) (M_PI*M_PI/(alpha*alpha));

229
230
    int firstz = (start == 0 ? 1 : 0);
    for (int kx = start; kx < end; kx++) {
231
        int mx = (kx < (gridx+1)/2) ? kx : kx-gridx;
peastman's avatar
peastman committed
232
        float mhx = mx*(float)recipBoxVectors[0][0];
233
234
235
        float bx = scaleFactor*bsplineModuli[0][kx];
        for (int ky = 0; ky < gridy; ky++) {
            int my = (ky < (gridy+1)/2) ? ky : ky-gridy;
peastman's avatar
peastman committed
236
            float mhy = mx*(float)recipBoxVectors[1][0] + my*(float)recipBoxVectors[1][1];
237
238
            float mhx2y2 = mhx*mhx + mhy*mhy;
            float bxby = bx*bsplineModuli[1][ky];
239
240
            for (int kz = firstz; kz < zsize; kz++) {
                int index = kx*yzsize + ky*zsize + kz;
241
                int mz = (kz < (gridz+1)/2) ? kz : kz-gridz;
peastman's avatar
peastman committed
242
                float mhz = mx*(float)recipBoxVectors[2][0] + my*(float)recipBoxVectors[2][1] + mz*(float)recipBoxVectors[2][2];
243
                float bz = bsplineModuli[2][kz];
244
245
                float m2 = mhx2y2 + mhz*mhz;
                float denom = m2*bxby*bz;
246
247
248
249
250
251
252
                recipEterm[index] = exp(-recipExpFactor*m2)/denom;
            }
            firstz = 0;
        }
    }
}

Peter Eastman's avatar
Peter Eastman committed
253
static double reciprocalEnergy(int start, int end, vector<complex<float> >& grid, vector<float>& recipEterm, int gridx, int gridy, int gridz, double alpha, vector<float>* bsplineModuli, Vec3* periodicBoxVectors, Vec3* recipBoxVectors) {
254
255
    const unsigned int zsizeHalf = gridz/2+1;
    const unsigned int yzsizeHalf = gridy*zsizeHalf;
256

257
    double energy = 0.0;
258
259
260
261
262

    int firstz = (start == 0 ? 1 : 0);
    for (int kx = start; kx < end; kx++) {
        for (int ky = 0; ky < gridy; ky++) {
            for (int kz = firstz; kz < gridz; kz++) {
263
264
265
266
267
268
269
270
271
272
273
                int kx1, ky1, kz1;
                if (kz >= gridz/2+1) {
                    kx1 = (kx == 0 ? kx : gridx-kx);
                    ky1 = (ky == 0 ? ky : gridy-ky);
                    kz1 = gridz-kz;
                }
                else {
                    kx1 = kx;
                    ky1 = ky;
                    kz1 = kz;
                }
274
                int index = kx1*yzsizeHalf + ky1*zsizeHalf + kz1;
Peter Eastman's avatar
Peter Eastman committed
275
276
                float gridReal = grid[index].real();
                float gridImag = grid[index].imag();
277
                energy += recipEterm[index]*(gridReal*gridReal+gridImag*gridImag);
278
279
280
281
            }
            firstz = 0;
        }
    }
282
    return 0.5*energy;
283
284
}

285

Peter Eastman's avatar
Peter Eastman committed
286
static double reciprocalDispersionEnergy(int start, int end, vector<complex<float> >& grid, const vector<float>& recipEterm, int gridx, int gridy, int gridz, double alpha, vector<float>* bsplineModuli, Vec3* periodicBoxVectors, Vec3* recipBoxVectors) {
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
    const unsigned int zsizeHalf = gridz/2+1;
    const unsigned int yzsizeHalf = gridy*zsizeHalf;

    double energy = 0.0;

    for (int kx = start; kx < end; kx++) {
        for (int ky = 0; ky < gridy; ky++) {
            for (int kz = 0; kz < gridz; kz++) {
                int kx1, ky1, kz1;
                if (kz >= gridz/2+1) {
                    kx1 = (kx == 0 ? kx : gridx-kx);
                    ky1 = (ky == 0 ? ky : gridy-ky);
                    kz1 = gridz-kz;
                }
                else {
                    kx1 = kx;
                    ky1 = ky;
                    kz1 = kz;
                }
                int index = kx1*yzsizeHalf + ky1*zsizeHalf + kz1;
Peter Eastman's avatar
Peter Eastman committed
307
308
                float gridReal = grid[index].real();
                float gridImag = grid[index].imag();
309
                energy += recipEterm[index]*(gridReal*gridReal+gridImag*gridImag);
310
311
312
            }
        }
    }
313
    return 0.5f*energy;
314
315
}

316

Peter Eastman's avatar
Peter Eastman committed
317
static void reciprocalConvolution(int start, int end, vector<complex<float> >& grid, vector<float>& recipEterm) {
318
319
    for (int index = start; index < end; index++) {
        float eterm = recipEterm[index];
Peter Eastman's avatar
Peter Eastman committed
320
        grid[index] *= eterm;
321
322
323
    }
}

Peter Eastman's avatar
Peter Eastman committed
324
static void interpolateForces(float* posq, vector<float>& force, vector<float>& grid, int gridx, int gridy, int gridz, int numParticles, Vec3* periodicBoxVectors, Vec3* recipBoxVectors, atomic<int>& atomicCounter, const float epsilonFactor, int numThreads) {
peastman's avatar
peastman committed
325
326
327
328
329
    fvec4 boxSize((float) periodicBoxVectors[0][0], (float) periodicBoxVectors[1][1], (float) periodicBoxVectors[2][2], 0);
    fvec4 invBoxSize((float) recipBoxVectors[0][0], (float) recipBoxVectors[1][1], (float) recipBoxVectors[2][2], 0);
    fvec4 recipBoxVec0((float) recipBoxVectors[0][0], (float) recipBoxVectors[0][1], (float) recipBoxVectors[0][2], 0);
    fvec4 recipBoxVec1((float) recipBoxVectors[1][0], (float) recipBoxVectors[1][1], (float) recipBoxVectors[1][2], 0);
    fvec4 recipBoxVec2((float) recipBoxVectors[2][0], (float) recipBoxVectors[2][1], (float) recipBoxVectors[2][2], 0);
330
331
332
333
    fvec4 gridSize(gridx, gridy, gridz, 0);
    ivec4 gridSizeInt(gridx, gridy, gridz, 0);
    fvec4 one(1);
    fvec4 scale(1.0f/(PME_ORDER-1));
334
335

    const int groupSize = max(1, numParticles / (10 * numThreads));
336
    while (true) {
337
338
        int start = atomicCounter.fetch_add(groupSize);
        if (start >= numParticles)
339
340
            break;

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
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
        int end = min(start + groupSize, numParticles);

        for (int i = start; i < end; i++) {
            // Find the position relative to the nearest grid point.

            fvec4 pos(&posq[4*i]);
            float posInBox[4];
            (pos-boxSize*floor(pos*invBoxSize)).store(posInBox);
            fvec4 t = posInBox[0]*recipBoxVec0 + posInBox[1]*recipBoxVec1 + posInBox[2]*recipBoxVec2;
            t = (t-floor(t))*gridSize;
            ivec4 ti = t;
            fvec4 dr = t-ti;
            ivec4 gridIndex = ti-(gridSizeInt&ti==gridSizeInt);

            // Compute the B-spline coefficients.

            fvec4 data[PME_ORDER];
            fvec4 ddata[PME_ORDER];
            data[PME_ORDER-1] = 0.0f;
            data[1] = dr;
            data[0] = one-dr;
            for (int j = 3; j < PME_ORDER; j++) {
                fvec4 div(1.0f/(j-1));
                data[j-1] = div*dr*data[j-2];
                for (int k = 1; k < j-1; k++)
                    data[j-k-1] = div*((dr+k)*data[j-k-2]+(fvec4(j-k)-dr)*data[j-k-1]);
                data[0] = div*(one-dr)*data[0];
            }
            ddata[0] = -data[0];
            for (int j = 1; j < PME_ORDER; j++)
                ddata[j] = data[j-1]-data[j];
            data[PME_ORDER-1] = scale*dr*data[PME_ORDER-2];
            for (int j = 1; j < (PME_ORDER-1); j++)
                data[PME_ORDER-j-1] = scale*((dr+j)*data[PME_ORDER-j-2]+(fvec4(PME_ORDER-j)-dr)*data[PME_ORDER-j-1]);
            data[0] = scale*(one-dr)*data[0];

            // Compute the force on this atom.

            int gridIndexX = gridIndex[0];
            int gridIndexY = gridIndex[1];
            int gridIndexZ = gridIndex[2];
            if (gridIndexX < 0)
                return; // This happens when a simulation blows up and coordinates become NaN.
            int zindex[PME_ORDER];
            for (int j = 0; j < PME_ORDER; j++) {
                zindex[j] = gridIndexZ+j;
                zindex[j] -= (zindex[j] >= gridz ? gridz : 0);
            }
            fvec4 zdata[PME_ORDER];
            for (int j = 0; j < PME_ORDER; j++)
                zdata[j] = fvec4(data[j][2], data[j][2], ddata[j][2], 0);
            fvec4 f = 0.0f;
            for (int ix = 0; ix < PME_ORDER; ix++) {
                int xbase = gridIndexX+ix;
                xbase -= (xbase >= gridx ? gridx : 0);
                xbase = xbase*gridy*gridz;
                float dx = data[ix][0];
                float ddx = ddata[ix][0];
                fvec4 xdata(ddx, dx, dx, 0);

                for (int iy = 0; iy < PME_ORDER; iy++) {
                    int ybase = gridIndexY+iy;
                    ybase -= (ybase >= gridy ? gridy : 0);
                    ybase = xbase + ybase*gridz;
                    float dy = data[iy][1];
                    float ddy = ddata[iy][1];
                    fvec4 xydata = xdata*fvec4(dy, ddy, dy, 0);

                    for (int iz = 0; iz < PME_ORDER; iz++) {
                        fvec4 gridValue(grid[ybase+zindex[iz]]);
                        f = f+xydata*zdata[iz]*gridValue;
                    }
413
414
                }
            }
415
416
417
418
419
420
            f *= -epsilonFactor*posq[4*i+3];
            float fc[4];
            f.store(fc);
            force[4*i+0] = fc[0]*gridx*(float)recipBoxVectors[0][0];
            force[4*i+1] = fc[0]*gridx*(float)recipBoxVectors[1][0]+fc[1]*gridy*(float)recipBoxVectors[1][1];
            force[4*i+2] = fc[0]*gridx*(float)recipBoxVectors[2][0]+fc[1]*gridy*(float)recipBoxVectors[2][1]+fc[2]*gridz*(float)recipBoxVectors[2][2];
421
422
423
424
        }
    }
}

425
static void* threadBody(void* args) {
426
427
    CpuCalcPmeReciprocalForceKernel& owner = *reinterpret_cast<CpuCalcPmeReciprocalForceKernel*>(args);
    owner.runMainThread();
428
429
430
    return 0;
}

431
void CpuCalcPmeReciprocalForceKernel::initialize(int xsize, int ysize, int zsize, int numParticles, double alpha, bool deterministic) {
432
433
    if (!hasInitializedThreads) {
        numThreads = getNumProcessors();
434
435
436
        char* threadsEnv = getenv("OPENMM_CPU_THREADS");
        if (threadsEnv != NULL)
            stringstream(threadsEnv) >> numThreads;
437
438
        hasInitializedThreads = true;
    }
439
    threadEnergy.resize(numThreads);
Peter Eastman's avatar
Peter Eastman committed
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
    gridx = findFFTDimension(xsize);
    gridy = findFFTDimension(ysize);
    gridz = findFFTDimension(zsize);
    gridShape.push_back(gridx);
    gridShape.push_back(gridy);
    gridShape.push_back(gridz);
    fftAxes.push_back(0);
    fftAxes.push_back(1);
    fftAxes.push_back(2);
    realGridStride.push_back(gridy*gridz*sizeof(float));
    realGridStride.push_back(gridz*sizeof(float));
    realGridStride.push_back(sizeof(float));
    complexGridStride.push_back(gridy*(gridz/2+1)*sizeof(complex<float>));
    complexGridStride.push_back((gridz/2+1)*sizeof(complex<float>));
    complexGridStride.push_back(sizeof(complex<float>));
455
456
    this->numParticles = numParticles;
    this->alpha = alpha;
457
    this->deterministic = deterministic;
458
    force.resize(4*numParticles);
459
    recipEterm.resize(gridx*gridy*gridz);
460
    
461
462
    // Initialize threads.
    
463
    isFinished = false;
464
    mainThread = thread(threadBody, this);
465
    
466
467
    // Wait until the main thread is up and running.
    
468
469
470
471
472
    {
        unique_lock<mutex> ul(lock);
        while (!isFinished)
            endCondition.wait(ul);
    }
473
    
Peter Eastman's avatar
Peter Eastman committed
474
475
476
477
    // Initialize the FFT grids.

    realGrids.resize(numThreads, vector<float>(gridx*gridy*gridz+3));
    complexGrid.resize(gridx*gridy*(gridz/2+1));
478
479
480
    
    // Initialize the b-spline moduli.

481
    int maxSize = std::max(std::max(gridx, gridy), gridz);
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
    vector<double> data(PME_ORDER);
    vector<double> ddata(PME_ORDER);
    vector<double> bsplinesData(maxSize);
    data[PME_ORDER-1] = 0.0;
    data[1] = 0.0;
    data[0] = 1.0;
    for (int i = 3; i < PME_ORDER; i++) {
        double div = 1.0/(i-1.0);
        data[i-1] = 0.0;
        for (int j = 1; j < (i-1); j++)
            data[i-j-1] = div*(j*data[i-j-2]+(i-j)*data[i-j-1]);
        data[0] = div*data[0];
    }

    // Differentiate.

    ddata[0] = -data[0];
    for (int i = 1; i < PME_ORDER; i++)
        ddata[i] = data[i-1]-data[i];
    double div = 1.0/(PME_ORDER-1);
    data[PME_ORDER-1] = 0.0;
    for (int i = 1; i < (PME_ORDER-1); i++)
        data[PME_ORDER-i-1] = div*(i*data[PME_ORDER-i-2]+(PME_ORDER-i)*data[PME_ORDER-i-1]);
    data[0] = div*data[0];
    for (int i = 0; i < maxSize; i++)
        bsplinesData[i] = 0.0;
    for (int i = 1; i <= PME_ORDER; i++)
        bsplinesData[i] = data[i-1];

    // Evaluate the actual bspline moduli for X/Y/Z.

    bsplineModuli[0].resize(gridx);
    bsplineModuli[1].resize(gridy);
    bsplineModuli[2].resize(gridz);
    for (int dim = 0; dim < 3; dim++) {
        int ndata = bsplineModuli[dim].size();
        vector<float>& moduli = bsplineModuli[dim];
        for (int i = 0; i < ndata; i++) {
            double sc = 0.0;
            double ss = 0.0;
            for (int j = 0; j < ndata; j++) {
                double arg = (2.0*M_PI*i*j)/ndata;
                sc += bsplinesData[j]*cos(arg);
                ss += bsplinesData[j]*sin(arg);
            }
            moduli[i] = (float) (sc*sc+ss*ss);
        }
        for (int i = 0; i < ndata; i++)
            if (moduli[i] < 1.0e-7f)
531
                moduli[i] = (moduli[(i-1+ndata)%ndata]+moduli[(i+1)%ndata])*0.5f;
532
533
534
    }
}

535
CpuCalcPmeReciprocalForceKernel::~CpuCalcPmeReciprocalForceKernel() {
536
    isDeleted = true;
537
538
539
540
    lock.lock();
    startCondition.notify_all();
    lock.unlock();
    mainThread.join();
541
542
}

543
544
void CpuCalcPmeReciprocalForceKernel::runMainThread() {
    // This is the main thread that coordinates all the other ones.
545

546
    unique_lock<mutex> ul(lock);
547
    isFinished = true;
548
    endCondition.notify_one();
549
550
551
552
    ThreadPool threads(numThreads);
    while (true) {
        // Wait for the signal to start.

553
        startCondition.wait(ul);
554
555
556
        if (isDeleted)
            break;
        posq = io->getPosq();
peastman's avatar
peastman committed
557
        atomicCounter = 0;
peastman's avatar
peastman committed
558
        threads.execute([&] (ThreadPool& threads, int threadIndex) { runWorkerThread(threads, threadIndex); }); // Signal threads to perform charge spreading.
559
560
561
        threads.waitForThreads();
        threads.resumeThreads(); // Signal threads to sum the charge grids.
        threads.waitForThreads();
Peter Eastman's avatar
Peter Eastman committed
562
        pocketfft::r2c(gridShape, realGridStride, complexGridStride, fftAxes, true, realGrids[0].data(), complexGrid.data(), 1.0f, 0);
563
564
565
566
567
568
569
        if (lastBoxVectors[0] != periodicBoxVectors[0] || lastBoxVectors[1] != periodicBoxVectors[1] || lastBoxVectors[2] != periodicBoxVectors[2]) {
            threads.resumeThreads(); // Signal threads to compute the reciprocal scale factors.
            threads.waitForThreads();
        }
        if (includeEnergy) {
            threads.resumeThreads(); // Signal threads to compute energy.
            threads.waitForThreads();
peastman's avatar
peastman committed
570
571
            for (auto e : threadEnergy)
                energy += e;
572
573
574
        }
        threads.resumeThreads(); // Signal threads to perform reciprocal convolution.
        threads.waitForThreads();
Peter Eastman's avatar
Peter Eastman committed
575
        pocketfft::c2r(gridShape, complexGridStride, realGridStride, fftAxes, false, complexGrid.data(), realGrids[0].data(), 1.0f, 0);
peastman's avatar
peastman committed
576
        atomicCounter = 0;
577
578
579
580
581
582
        threads.resumeThreads(); // Signal threads to interpolate forces.
        threads.waitForThreads();
        isFinished = true;
        lastBoxVectors[0] = periodicBoxVectors[0];
        lastBoxVectors[1] = periodicBoxVectors[1];
        lastBoxVectors[2] = periodicBoxVectors[2];
583
        endCondition.notify_one();
584
    }
585
586
}

587
588
589
590
591
592
void CpuCalcPmeReciprocalForceKernel::runWorkerThread(ThreadPool& threads, int index) {
    int gridxStart = (index*gridx)/numThreads;
    int gridxEnd = ((index+1)*gridx)/numThreads;
    int gridSize = (gridx*gridy*gridz+3)/4;
    int gridStart = 4*((index*gridSize)/numThreads);
    int gridEnd = 4*(((index+1)*gridSize)/numThreads);
593
    int complexSize = gridx*gridy*(gridz/2+1);
594
    int complexStart = std::max(1, ((index*complexSize)/numThreads));
595
    int complexEnd = (((index+1)*complexSize)/numThreads);
596
    const float epsilonFactor = sqrt(ONE_4PI_EPS0);
Peter Eastman's avatar
Peter Eastman committed
597
    spreadCharge(posq, realGrids[index], gridx, gridy, gridz, numParticles, periodicBoxVectors, recipBoxVectors, atomicCounter, epsilonFactor, index, numThreads, deterministic);
598
    threads.syncThreads();
Peter Eastman's avatar
Peter Eastman committed
599
    int numGrids = realGrids.size();
600
    for (int i = gridStart; i < gridEnd; i += 4) {
Peter Eastman's avatar
Peter Eastman committed
601
        fvec4 sum(&realGrids[0][i]);
602
        for (int j = 1; j < numGrids; j++)
Peter Eastman's avatar
Peter Eastman committed
603
604
            sum += fvec4(&realGrids[j][i]);
        sum.store(&realGrids[0][i]);
605
    }
606
    threads.syncThreads();
607
608
    if (lastBoxVectors[0] != periodicBoxVectors[0] || lastBoxVectors[1] != periodicBoxVectors[1] || lastBoxVectors[2] != periodicBoxVectors[2]) {
        computeReciprocalEterm(gridxStart, gridxEnd, gridx, gridy, gridz, recipEterm, alpha, bsplineModuli, periodicBoxVectors, recipBoxVectors);
609
        threads.syncThreads();
610
611
612
    }
    if (includeEnergy) {
        threadEnergy[index] = reciprocalEnergy(gridxStart, gridxEnd, complexGrid, recipEterm, gridx, gridy, gridz, alpha, bsplineModuli, periodicBoxVectors, recipBoxVectors);
613
614
        threads.syncThreads();
    }
615
616
    reciprocalConvolution(complexStart, complexEnd, complexGrid, recipEterm);
    threads.syncThreads();
Peter Eastman's avatar
Peter Eastman committed
617
    interpolateForces(posq, force, realGrids[0], gridx, gridy, gridz, numParticles, periodicBoxVectors, recipBoxVectors, atomicCounter, epsilonFactor, numThreads);
618
619
}

peastman's avatar
peastman committed
620
void CpuCalcPmeReciprocalForceKernel::beginComputation(IO& io, const Vec3* periodicBoxVectors, bool includeEnergy) {
621
    this->io = &io;
peastman's avatar
peastman committed
622
623
624
    this->periodicBoxVectors[0] = periodicBoxVectors[0];
    this->periodicBoxVectors[1] = periodicBoxVectors[1];
    this->periodicBoxVectors[2] = periodicBoxVectors[2];
625
626
    this->includeEnergy = includeEnergy;
    energy = 0.0;
627
    ReferenceForce::invertBoxVectors(periodicBoxVectors, recipBoxVectors);
peastman's avatar
peastman committed
628
629
630

    // Do the calculation.

631
    unique_lock<mutex> ul(lock);
632
    isFinished = false;
633
    startCondition.notify_one();
634
635
}

636
double CpuCalcPmeReciprocalForceKernel::finishComputation(IO& io) {
637
638
639
640
641
    {
        unique_lock<mutex> ul(lock);
        while (!isFinished) {
            endCondition.wait(ul);
        }
642
643
    }
    io.setForce(&force[0]);
644
    return energy;
peastman's avatar
peastman committed
645
}
646
647

bool CpuCalcPmeReciprocalForceKernel::isProcessorSupported() {
648
    return isVec4Supported();
649
}
650

651
652
653
654
655
656
657
void CpuCalcPmeReciprocalForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
    alpha = this->alpha;
    nx = gridx;
    ny = gridy;
    nz = gridz;
}

Peter Eastman's avatar
Peter Eastman committed
658
int CpuCalcPmeReciprocalForceKernel::findFFTDimension(int minimum) {
659
660
661
662
663
664
    if (minimum < 1)
        return 1;
    while (true) {
        // Attempt to factor the current value.

        int unfactored = minimum;
Peter Eastman's avatar
Peter Eastman committed
665
        for (int factor = 2; factor < 9; factor++) {
666
667
668
            while (unfactored > 1 && unfactored%factor == 0)
                unfactored /= factor;
        }
Peter Eastman's avatar
Peter Eastman committed
669
        if (unfactored == 1 || unfactored == 11)
670
671
672
673
            return minimum;
        minimum++;
    }
}
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699

/*
 * Everything below here is just a clone of the above, but to handle the dispersion term
 * instead of electrostatics.
 */

bool CpuCalcPmeReciprocalForceKernel::hasInitializedThreads = false;
int CpuCalcPmeReciprocalForceKernel::numThreads = 0;


class CpuCalcDispersionPmeReciprocalForceKernel::ComputeTask : public ThreadPool::Task {
public:
    ComputeTask(CpuCalcDispersionPmeReciprocalForceKernel& owner) : owner(owner) {
    }
    void execute(ThreadPool& threads, int threadIndex) {
        owner.runWorkerThread(threads, threadIndex);
    }
    CpuCalcDispersionPmeReciprocalForceKernel& owner;
};

static void* dispersionThreadBody(void* args) {
    CpuCalcDispersionPmeReciprocalForceKernel& owner = *reinterpret_cast<CpuCalcDispersionPmeReciprocalForceKernel*>(args);
    owner.runMainThread();
    return 0;
}

700
void CpuCalcDispersionPmeReciprocalForceKernel::initialize(int xsize, int ysize, int zsize, int numParticles, double alpha, bool deterministic) {
701
702
703
704
705
706
707
708
    if (!hasInitializedThreads) {
        numThreads = getNumProcessors();
        char* threadsEnv = getenv("OPENMM_CPU_THREADS");
        if (threadsEnv != NULL)
            stringstream(threadsEnv) >> numThreads;
        hasInitializedThreads = true;
    }
    threadEnergy.resize(numThreads);
Peter Eastman's avatar
Peter Eastman committed
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
    gridx = findFFTDimension(xsize);
    gridy = findFFTDimension(ysize);
    gridz = findFFTDimension(zsize);
    gridShape.push_back(gridx);
    gridShape.push_back(gridy);
    gridShape.push_back(gridz);
    fftAxes.push_back(0);
    fftAxes.push_back(1);
    fftAxes.push_back(2);
    realGridStride.push_back(gridy*gridz*sizeof(float));
    realGridStride.push_back(gridz*sizeof(float));
    realGridStride.push_back(sizeof(float));
    complexGridStride.push_back(gridy*(gridz/2+1)*sizeof(complex<float>));
    complexGridStride.push_back((gridz/2+1)*sizeof(complex<float>));
    complexGridStride.push_back(sizeof(complex<float>));
724
725
    this->numParticles = numParticles;
    this->alpha = alpha;
726
    this->deterministic = deterministic;
727
728
729
730
731
732
    force.resize(4*numParticles);
    recipEterm.resize(gridx*gridy*gridz);
    
    // Initialize threads.
    
    isFinished = false;
733
    mainThread = thread(dispersionThreadBody, this);
734
735
736
    
    // Wait until the main thread is up and running.
    
737
738
739
740
741
    {
        unique_lock<mutex> ul(lock);
        while (!isFinished)
            endCondition.wait(ul);
    }
Peter Eastman's avatar
Peter Eastman committed
742
743
744
745
746

    // Initialize the FFT grids.

    realGrids.resize(numThreads, vector<float>(gridx*gridy*gridz+3));
    complexGrid.resize(gridx*gridy*(gridz/2+1));
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
    
    // Initialize the b-spline moduli.

    int maxSize = std::max(std::max(gridx, gridy), gridz);
    vector<double> data(PME_ORDER);
    vector<double> ddata(PME_ORDER);
    vector<double> bsplinesData(maxSize);
    data[PME_ORDER-1] = 0.0;
    data[1] = 0.0;
    data[0] = 1.0;
    for (int i = 3; i < PME_ORDER; i++) {
        double div = 1.0/(i-1.0);
        data[i-1] = 0.0;
        for (int j = 1; j < (i-1); j++)
            data[i-j-1] = div*(j*data[i-j-2]+(i-j)*data[i-j-1]);
        data[0] = div*data[0];
    }

    // Differentiate.

    ddata[0] = -data[0];
    for (int i = 1; i < PME_ORDER; i++)
        ddata[i] = data[i-1]-data[i];
    double div = 1.0/(PME_ORDER-1);
    data[PME_ORDER-1] = 0.0;
    for (int i = 1; i < (PME_ORDER-1); i++)
        data[PME_ORDER-i-1] = div*(i*data[PME_ORDER-i-2]+(PME_ORDER-i)*data[PME_ORDER-i-1]);
    data[0] = div*data[0];
    for (int i = 0; i < maxSize; i++)
        bsplinesData[i] = 0.0;
    for (int i = 1; i <= PME_ORDER; i++)
        bsplinesData[i] = data[i-1];

    // Evaluate the actual bspline moduli for X/Y/Z.

    bsplineModuli[0].resize(gridx);
    bsplineModuli[1].resize(gridy);
    bsplineModuli[2].resize(gridz);
    for (int dim = 0; dim < 3; dim++) {
        int ndata = bsplineModuli[dim].size();
        vector<float>& moduli = bsplineModuli[dim];
        for (int i = 0; i < ndata; i++) {
            double sc = 0.0;
            double ss = 0.0;
            for (int j = 0; j < ndata; j++) {
                double arg = (2.0*M_PI*i*j)/ndata;
                sc += bsplinesData[j]*cos(arg);
                ss += bsplinesData[j]*sin(arg);
            }
            moduli[i] = (float) (sc*sc+ss*ss);
        }
        for (int i = 0; i < ndata; i++)
            if (moduli[i] < 1.0e-7f)
                moduli[i] = (moduli[i-1]+moduli[i+1])*0.5f;
    }
}

CpuCalcDispersionPmeReciprocalForceKernel::~CpuCalcDispersionPmeReciprocalForceKernel() {
    isDeleted = true;
806
807
808
809
    lock.lock();
    startCondition.notify_all();
    lock.unlock();
    mainThread.join();
810
811
812
813
814
}

void CpuCalcDispersionPmeReciprocalForceKernel::runMainThread() {
    // This is the main thread that coordinates all the other ones.

815
    unique_lock<mutex> ul(lock);
816
    isFinished = true;
817
    endCondition.notify_one();
818
819
820
821
    ThreadPool threads(numThreads);
    while (true) {
        // Wait for the signal to start.

822
        startCondition.wait(ul);
823
824
825
826
        if (isDeleted)
            break;
        posq = io->getPosq();
        ComputeTask task(*this);
peastman's avatar
peastman committed
827
        atomicCounter = 0;
828
829
830
831
        threads.execute(task); // Signal threads to perform charge spreading.
        threads.waitForThreads();
        threads.resumeThreads(); // Signal threads to sum the charge grids.
        threads.waitForThreads();
Peter Eastman's avatar
Peter Eastman committed
832
        pocketfft::r2c(gridShape, realGridStride, complexGridStride, fftAxes, true, realGrids[0].data(), complexGrid.data(), 1.0f, 0);
833
834
835
836
837
838
839
        if (lastBoxVectors[0] != periodicBoxVectors[0] || lastBoxVectors[1] != periodicBoxVectors[1] || lastBoxVectors[2] != periodicBoxVectors[2]) {
            threads.resumeThreads(); // Signal threads to compute the reciprocal scale factors.
            threads.waitForThreads();
        }
        if (includeEnergy) {
            threads.resumeThreads(); // Signal threads to compute energy.
            threads.waitForThreads();
peastman's avatar
peastman committed
840
841
            for (auto e : threadEnergy)
                energy += e;
842
843
844
        }
        threads.resumeThreads(); // Signal threads to perform reciprocal convolution.
        threads.waitForThreads();
Peter Eastman's avatar
Peter Eastman committed
845
        pocketfft::c2r(gridShape, complexGridStride, realGridStride, fftAxes, false, complexGrid.data(), realGrids[0].data(), 1.0f, 0);
peastman's avatar
peastman committed
846
        atomicCounter = 0;
847
848
849
850
851
852
        threads.resumeThreads(); // Signal threads to interpolate forces.
        threads.waitForThreads();
        isFinished = true;
        lastBoxVectors[0] = periodicBoxVectors[0];
        lastBoxVectors[1] = periodicBoxVectors[1];
        lastBoxVectors[2] = periodicBoxVectors[2];
853
        endCondition.notify_one();
854
855
856
857
858
859
860
861
862
863
864
865
866
    }
}

void CpuCalcDispersionPmeReciprocalForceKernel::runWorkerThread(ThreadPool& threads, int index) {
    int gridxStart = (index*gridx)/numThreads;
    int gridxEnd = ((index+1)*gridx)/numThreads;
    int gridSize = (gridx*gridy*gridz+3)/4;
    int gridStart = 4*((index*gridSize)/numThreads);
    int gridEnd = 4*(((index+1)*gridSize)/numThreads);
    int complexSize = gridx*gridy*(gridz/2+1);
    int complexStart = std::max(1, ((index*complexSize)/numThreads));
    int complexEnd = (((index+1)*complexSize)/numThreads);
    const float epsilonFactor = 1.0f;
Peter Eastman's avatar
Peter Eastman committed
867
    spreadCharge(posq, realGrids[index], gridx, gridy, gridz, numParticles, periodicBoxVectors, recipBoxVectors, atomicCounter, epsilonFactor, index, numThreads, deterministic);
868
    threads.syncThreads();
Peter Eastman's avatar
Peter Eastman committed
869
    int numGrids = realGrids.size();
870
    for (int i = gridStart; i < gridEnd; i += 4) {
Peter Eastman's avatar
Peter Eastman committed
871
        fvec4 sum(&realGrids[0][i]);
872
        for (int j = 1; j < numGrids; j++)
Peter Eastman's avatar
Peter Eastman committed
873
874
            sum += fvec4(&realGrids[j][i]);
        sum.store(&realGrids[0][i]);
875
876
877
878
879
880
881
882
883
884
885
886
887
888
    }
    threads.syncThreads();
    if (lastBoxVectors[0] != periodicBoxVectors[0] || lastBoxVectors[1] != periodicBoxVectors[1] || lastBoxVectors[2] != periodicBoxVectors[2]) {
        computeReciprocalDispersionEterm(gridxStart, gridxEnd, gridx, gridy, gridz, recipEterm, alpha, bsplineModuli, periodicBoxVectors, recipBoxVectors);
        threads.syncThreads();
    }
    if (includeEnergy) {
        threadEnergy[index] = reciprocalDispersionEnergy(gridxStart, gridxEnd, complexGrid, recipEterm, gridx, gridy, gridz, alpha, bsplineModuli, periodicBoxVectors, recipBoxVectors);
        threads.syncThreads();
    }
    // For dispersion, we include the {0,0,0} term, so the start point needs to be redefined
    complexStart = (index*complexSize)/numThreads;
    reciprocalConvolution(complexStart, complexEnd, complexGrid, recipEterm);
    threads.syncThreads();
Peter Eastman's avatar
Peter Eastman committed
889
    interpolateForces(posq, force, realGrids[0], gridx, gridy, gridz, numParticles, periodicBoxVectors, recipBoxVectors, atomicCounter, epsilonFactor, numThreads);
890
891
}

892
void CpuCalcDispersionPmeReciprocalForceKernel::beginComputation(CalcPmeReciprocalForceKernel::IO& io, const Vec3* periodicBoxVectors, bool includeEnergy) {
893
894
895
896
897
898
    this->io = &io;
    this->periodicBoxVectors[0] = periodicBoxVectors[0];
    this->periodicBoxVectors[1] = periodicBoxVectors[1];
    this->periodicBoxVectors[2] = periodicBoxVectors[2];
    this->includeEnergy = includeEnergy;
    energy = 0.0;
899
    ReferenceForce::invertBoxVectors(periodicBoxVectors, recipBoxVectors);
900
901
902

    // Do the calculation.

903
    unique_lock<mutex> ul(lock);
904
    isFinished = false;
905
    startCondition.notify_one();
906
907
}

908
double CpuCalcDispersionPmeReciprocalForceKernel::finishComputation(CalcPmeReciprocalForceKernel::IO& io) {
909
910
911
912
913
    {
        unique_lock<mutex> ul(lock);
        while (!isFinished) {
            endCondition.wait(ul);
        }
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
    }
    io.setForce(&force[0]);
    return energy;
}

bool CpuCalcDispersionPmeReciprocalForceKernel::isProcessorSupported() {
    return isVec4Supported();
}

void CpuCalcDispersionPmeReciprocalForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
    alpha = this->alpha;
    nx = gridx;
    ny = gridy;
    nz = gridz;
}

Peter Eastman's avatar
Peter Eastman committed
930
int CpuCalcDispersionPmeReciprocalForceKernel::findFFTDimension(int minimum) {
931
932
933
934
935
936
    if (minimum < 1)
        return 1;
    while (true) {
        // Attempt to factor the current value.

        int unfactored = minimum;
Peter Eastman's avatar
Peter Eastman committed
937
        for (int factor = 2; factor < 9; factor++) {
938
939
940
            while (unfactored > 1 && unfactored%factor == 0)
                unfactored /= factor;
        }
Peter Eastman's avatar
Peter Eastman committed
941
        if (unfactored == 1 || unfactored == 11)
942
943
944
945
            return minimum;
        minimum++;
    }
}