Interacting Particles

In this section we are going to take a shot at the problem of interacting particles. A typical application of this problem is to the motion of charged plasma particles that interact electrostatically with each other. Plasmas are actually even more complex, because accelerated particles radiate and some of this radiation may be captured by other particles, and this adds enormous complexity to the models. Other applications include molecular dynamics and, right at the other end of the spectrum, galactic dynamics.

As far as MPI is concerned, we are going to take another look at defining MPI types and at collective communications. We are going to arrange all particles into an array and the way that this array is going to be partitioned between processes of the computational pool is very reminiscent of how MPI-IO partitions files. So this section may be thought of as a preparation to MPI-IO as well.

Let me begin by reminding you about the Coulomb law. If you have two electrically charged particles, with charges q₁ and q₂, then the force with which particle 2 acts on particle 1 is given by

 

$$- \frac{1}{4\pi \epsilon_0} \frac{q_1 q_2}{r_{12}^2} \frac{\boldsymbol{r}_{12}}{r_{12}}\quad\textrm{where}$$ F₁₂ = (5.4)

$$\boldsymbol{r}_2 - \boldsymbol{r}_1\quad\textrm{and}$$ r₁₂ =

$$\sqrt{\left(r^x_{12}\right)^2 + \left(r^y_{12}\right)^2 + \left(r^z_{12}\right)^2}$$ r₁₂ =

where $\epsilon_0$ is the dialectric constant of vacuum and vector r₁₂ points from particle 1 to particle 2 (this is easy to see - simply assume that particle 1 sits in the middle of the world so that r₁ = 0). If the particle charges are of the opposite sign, then q₁ q₂ is negative. This cancels the minus in front of $1/(4\pi \epsilon_0)$ and particle 1 ends up being accelerated towards particle 2. The particles attract each other. If the particles have charges of the same sign, q₁ q₂ > 0 and then particle 1 is accelerated in the direction that is opposite to where particle 2 is. The particles repel each other.

The force with which particle 1 acts on particle 2 is the opposite to F₁₂, i.e.,

F₂₁ = - F₁₂

Sometimes equation (5.4) is printed without the minus in front of $1/(4\pi \epsilon_0)$. Whether the minus should or should not be there depends on the definitions of F₁₂ and r₁₂. If r₁₂ is the vector that points from particle 1 to particle 2, i.e., r₂ - r₁, and if F₁₂ is the force with which particle 2 acts on particle 1, then the minus should be there. Physics is one of these somewhat troubling disciplines, where you need to understand what happens and where you need to be very precise about what is what, otherwise it's easy to get confused.

Consider N particles of various charges scattered at random in some region of space. In order to evaluate the total force that acts on particle number i we have to evaluate the sum:

$$\boldsymbol{F}_i = - \frac{1}{4\pi\epsilon_0} \sum_{{j = 1... ... \frac{q_i q_j}{r_{ij}^2} \frac{\boldsymbol{r}_{ij}}{r_{ij}}$$ (5.5)

This sum has N-1 terms in it, because we assume that particle i does not act on itself. In order to evaluate forces that act on all N particles we have to evaluate N (N - 1) such terms. We say that this computation is going to scale like $\mathcal{O}(N^2)$. This is going to turn very expensive if N gets very large.

We are going to distribute the computation over n processes by making each process evaluate the force for N/n particles. Each process is therefore going to evaluate N/n (N - 1) terms, which may result in considerable speed up of the computation for a sufficiently large value of n. In particular, if we could make n = N, the computation would scale as $\mathcal{O}(N)$ instead of $\mathcal{O}(N^2)$. This is how nature goes about ``parallelizing'' this problem. You can think of each particle as being a separate processor that calculates its own evolution.

Within this example program we are going to push all N particles for 20 short time steps using a very simplistic numerical procedure that corresponds to the following equations of motion for particle i:

$$\begin{eqnarray*}m_i \frac{\mathrm{d} \boldsymbol{v}_i}{\mathrm{d} t} &=& - \fr... ...ac{\mathrm{d}\boldsymbol{r}_i}{\mathrm{d}t} &=& \boldsymbol{v}_i \end{eqnarray*}$$

We can always choose to use such time units in our computation that

$$\frac{1}{4\pi\epsilon_0} = 1$$

We are going to discretize the equations of motion by replacing dt with a final time interval $\Delta t$, so that:

$$\begin{eqnarray*}\boldsymbol{v}_i(t + \Delta t) &=& \boldsymbol{v}_i(t) - \frac... ...symbol{v}_i(t) + \boldsymbol{v}_i(t + \Delta t)} {2} \Delta t \end{eqnarray*}$$