# OpenMP: parallelizing calculation of pair interaction forces

I recently started to learn OpenMP. Albeit I have developed some intuition I still have some doubts on how to proceed under certain circumstances that are very useful for computational physicists.

My question: In molecular simulation is typical to find central forces, $$F=F(r)$$ where $$r = |\vec{r}_i -\vec{r}_j|$$ for particles $$i$$ and $$j$$. In addition, we can consider Newton's third law and then for a system of $$N$$ interacting particles we only need to run over half of the particles, since the other is simply obtained changing the sign. In C/C++$$^{(1)}$$:

 for(int i=0;i<N;i++){
for(int j=i+1;j<N;j++) {
F = force(r);
f[i] += F;
f[j] -= F;
}
}


This is a simple example where is necessary to sum and to subtract on a given vector and similar examples are very common even if no central forces are involved. However, I am not sure how to parallelize efficiently the aforementioned loops. One option could be to calculate everything instead of only one half:

#pragma omp for private(F), reduction(+:f[:N])
for(int i=0;i<N;i++){
for(int j=0;j<N;j++) {
if(i!=j){
F = force(x);
f[i] += F;
f[j] -= F;
}
}
}


But that is far from ideal since the calculation of F can be quite expensive and it will be done twice. The main problem that I see to parallelize the first code is the reduction since it'd have to be done for subtraction, reduction(-:f[:N]), and addition, reduction(+:f[:N]) with the ensuing conflict.

My alternative is through the use of two intermediate arrays, one for the sum component and other for the subtraction, lets call them sums_f and subs_f:

double *sums_f, *subs_f;
// Reserve memory for sums_f and subs_f and set them to zero
#pragma omp parallel
{
#pragma omp for reduction(+: sums_f[:N], -:subs_f[:N])
for(int i=0;i<N;i++){
for(int j=i+1;j<N;j++) {
F = force(r);
sums_f[i] += F;
subs_f[j] -= F;
}
}
}

#pragma omp for
for(int i=0;i<N;i++)
f[i] += sums_f[i] - subs_f[i];
}



But I think that can be very costly in terms of memory for large $$N$$.

In summary, I do no see an alternative to make them process efficient, without requiring several variables or additional steps that can downgrade the efficiency of parallelization. Do you have any suggestion on to proceed for such cases?

Notes: [1] For a more detailed explanation see: How to compute forces in multi-particle MD

20210513.1: IMPORTANT: Actually dynamically-allocated arrays are not suitable for reduction. This probably excludes its use in this context.

20210514.1: For object in 3D I use to store the force data per particle as a 3D dimensional array of $$N$$ rows by 3 columns. However, this means that is necessary to do 3x2 atomic steps. The result is that the gain in time of parallelization is lost. A better solution is to store the force data as a 1D array.

• Have you tried #pragma omp atomic for your addition and subtraction operations? – Charlie S May 12 at 0:27
• Why not do f[j] += -F; if you're worried about having to state the reduction twice? – Wolfgang Bangerth May 12 at 2:23
• I assume in the third code you means subsf[ j ] -=F; (note index, j not i), and the reduction also looks wrong – Ian Bush May 12 at 9:36
• By using Newton 3 you are saying that the force evaluations are no longer totally independent. As the operations are no longer totally independent there will almost certainly be an impact on parallelisation, as the base assumption of parallelisation is that the operations are independent. Thus you will almost certainly have to include some extra synchronisation or "dummy" variables. Personally as the increase in memory is linear (as opposed to quadratic or worse) I would not worry too much, and compare that with approaches that use atomic and similar to see what suits your needs best – Ian Bush May 12 at 9:41
• "Actually dynamically-allocated arrays are not suitable for reduction" - Really!!? It's certainly fine in Fortran. But you can always do the reduction "manually" - create a private array for each thread, sum into the private array, and then sum over the private arrays after the end of the main force loop. After all, that's probably what the OpenMP implementation does internally for a reduction. – Ian Bush May 14 at 19:20

So I decided to implement this, and of course I found a bunch of factors that I didn't foresaw in my earlier answer, which I deleted because it contained some patent nonsense.

Sequential code

      for (int ip=0; ip<N; ip++) {
for (int jp=ip+1; jp<N; jp++) {
struct force f = force_calc(points[ip],points[jp]);
sub_force( forces+jp,f );
}
}


The add/sub routines are to account for the different directions of the force arrow.

Parallel by looping over the full N^2 and ignore the extra work

    for (int step=0; step<STEPS; step++) {
for (int ip=0; ip<N; ip++) {
double sumx=0., sumy=0., sumf=0.;
#pragma omp parallel for reduction(+:sumx,sumy,sumf)
for (int jp=0; jp<N; jp++) {
if (ip==jp) continue;
struct force f = force_calc(points[ip],points[jp]);
sumx += f.x; sumy += f.y; sumf += f.f;
} // end parallel jp loop
struct force sumforce = {sumx,sumy,sumf};
} // end ip loop


Note that we can only parallelize the inner loop. If you try to parallelize both and add a collapse clause you will get incorrect results.

Parallel code 2, looping over the triangle, but using atomics in the sub_force routine:

#pragma omp parallel for schedule(guided,4)
for (int ip=0; ip<N; ip++) {
for (int jp=ip+1; jp<N; jp++) {
struct force f = force_calc(points[ip],points[jp]);
sub_force( forces+jp,f );
}
}


And now for the timings. I'm using a dual socket Cascade Lake with 56 cores total.

================ #threads = 1 ================
Sequential: 1.167970e+01; total force: 7.382722e+08
Full loop Parallel: 3.132034e+00; total force: 7.382722e+08, speedup= 3.73
Atomic update parallel: 1.136083e+01; total force: 7.382722e+08, speedup= 1.03
Sequential: 1.166757e+01; total force: 7.589678e+08
Full loop Parallel: 1.910710e+00; total force: 7.589678e+08, speedup= 6.11
Atomic update parallel: 1.733956e+00; total force: 7.588650e+08, speedup= 6.73
Sequential: 1.167301e+01; total force: 7.378324e+08
Full loop Parallel: 4.362817e+00; total force: 7.378324e+08, speedup= 2.68
Atomic update parallel: 1.087576e+00; total force: 7.376554e+08, speedup=10.73
Sequential: 1.168075e+01; total force: 7.428077e+08
Full loop Parallel: 9.035676e+00; total force: 7.428077e+08, speedup= 1.29
Atomic update parallel: 9.829427e-01; total force: 7.425870e+08, speedup=11.88


First big surprise: the double loop is 4x faster, (not 2x slower because of the double work), in the case without parallelism. That's because the triangular version jumps all through memory, meaning worse cache and TLB behavior.

Surprise two: the double loop gets worse and worse. I'm guessing that the array gets too short (10k particles) to get good efficiency.

Surprise three: using atomic updates is actually efficient. We're getting a decent speedup, although far from linear. I'm using a guided schedule because the work is not constant per ip iteration. Maybe if the force calculation is more expensive this will scale even better.

EDIT but for the best solution, we make everything atomic, and collapse the two loops. Now there is a lot more work to keep all the cores happy.

#pragma omp parallel for collapse(2)
for (int ip=0; ip<N; ip++) {
for (int jp=0; jp<N; jp++) {
if (ip==jp) continue;
struct force f = force_calc(points[ip],points[jp]);
} // end parallel jp loop
} // end ip loop


where the add routine using #omp atomic updates.

Timing:

================ #threads = 1 ================
Sequential: 2.029093e+01;
Full loop atomic: 3.114875e+01; speedup= 0.65
Full loop atomic: 1.739044e+00; speedup=11.67
Full loop atomic: 8.486451e-01; speedup=23.95
Full loop atomic: 5.630970e-01; speedup=36.09


This came as a bit of a surprise to me. Modern hardware must really know how to deal with that atomic stuff.

• Eijhout, Nice! thank you very much for you work, much appreciated. I have not internalized yet the schedule clause but I see is worthy to study it. I deduce that forces is an array of structs with components x,y,z. Is this approacht more efficient than using a 2D array in sequential? I ask because one of the nice things of OpenMP is that one does not have to big changes in the already-written codes. Since in parallel one can gain a factor 10x it might be interesting to design my codes as arrays of structs in the future. Thanks again. – CLP May 18 at 10:21
• 1. The schedule clause is needed here because the work per iteration differs, so you want some form of dynamic scheduling. Usually "guided" is cheaper than "dynamic". 2. Yes, I'm using an array of structures for cleanliness of coding, but I think there is no difference in efficiency. In fact, here it may be the best solution. 3. Raaahh!! Immediately after I posted this, I found the best solution: fully atomic, and collapse the two loops. See my edit. – Victor Eijkhout May 18 at 12:49

Not necessarily an answer, just easier to format code here. The #pragma omp atomic directive ensures that race conditions are avoided in a parallel loop. It is similar to critical but only works for simple expressions, like +=, and is faster. I believe this is the simplest way of achieving what you have in mind.

#pragma omp parallel for
for (int i = 0; i < N; i++)
{
for (int j = i + 1; j < N; j++)
{
F = force(r);
#pragma omp atomic
f[i] += F;
#pragma omp atomic
f[j] -= F;
}
}


Edit: now with more cleverness

If you want to spread the work of force(r) more evenly, define $$k = j(j+1)/2 - i - 1$$ $$j = \lceil\sqrt{(2k+2.25)}-0.5\rceil$$ $$i = j(j+1)/2 - k - 1$$

Don't ask me to prove that it creates a unique pair -- this was the result of messing around and trying stuff until it worked. I have tested it up to N=45000 (any more than you need long). This should avoid having to explicitly create the (i,j) pairs.

std::pair<int, int> ij(int k)
{
int j = std::ceil(sqrt(2 * k + 2.25) - 0.5);
int tri_j = j * (j + 1) / 2;
int i = tri_j - k - 1;
return std::make_pair(i, j);
}

//number of unique pairs (i,j)
int M = N*(N-1)/2;

#pragma omp parallel for
for (int k = 0; k < M; k++)
{
std::pair<int,int> ij = calc_ij(k);
int i = ij.first;
int j = ij.second;

//presumably the bottleneck
double F = force(r);

//practically instantaneous
#pragma omp atomic
f[i] += F;
#pragma omp atomic
f[j] -= F;
}



I'd be interested to see if the linear indexing improves execution time.

• Thanks for the code. As far as I know using atomic/critical inside the loop is a slower strategy than incorporating them outside the loop using a third variable, with both the loop and the final sumation inside a #pragma omp parallel (similar to third example)? Correct me if I am wrong. – CLP May 12 at 16:38
• Only timing will tell. That said, atomic is pretty fast and isn't the bottleneck (I assume force(r) consumes most of the computer's time). I would be more interested in combining i and j into a single loop so the work is more evenly distributed. That is: for(int k = 0; k < N*(N-1)/2; k++)... – Charlie S May 12 at 17:28
• Thanks for the update in the code, I did not know the std::pair tool, much appreciated. I wanted to do a trial, but I realized I had explored something similar: Note that since $M~N^2$, even for a "moderately large" number of particles the memory necessary is huge. In similar approaches I have tried, the computer freeze due to the little memory available – CLP May 13 at 22:39