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 <PERSON>'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. 603
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: https://scicomp.stackexchange.com/questions/28484/how-to-compute-forces-in-multi-particle-md
Updates:
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. 603