6.5 Particles and Fields
In many cases the interparticle potential may be safely neglected beyond
a few particle diameters. Defining a cutoff distance
$r_{co}$, we may estimate the importance of the neglected
"tail" of a potential $u(r)$ by
$
\int \limits_{\textstyle r_{co}}^{ \infty} u(r) 4 \pi r^{2} dr
$
If the decay of $u(r)$ is slow, this integral may be nonnegligible.
For example, the interaction between charged particles
decays only as $r^{1}$, and any cutoff will introduce a considerable
error. Another instance is astrophysics, where the gravitational potential
reaches too far to permit a simple cutoff.
For ionic systems that are globally neutral the Ewald summation
method known from solid state physics may be invoked. In the case of
hot plasmas or gravitating systems the ParticleMesh
methods are appropriate.
Subsections
6.5.1 Ewald summation
Let
$
u_{qq}=\frac{\textstyle q_{1}q_{2}}{\textstyle r}
$
be the ionion interaction between charged particles.
In the Ewald summation approach the basic cell containing
$N/2$ each of positive and negative charges in some
spatial arrangement is interpreted as a single crystallographic
element surrounded by an infinite number of identical copies of itself.
The entire system is then neutral and contains an infinite number of charges
situated at points $ r_{j+}$ and $r_{j}$, respectively. The
total potential at the position of some ion $i$ residing in the basic cell
is given by the finite difference of two infinite, diverging series:
$
\phi( r_{i})=
q \sum \limits_{j+}^{\infty}
\frac{\textstyle 1}{\textstyle  r_{i} r_{j+} } 
q \sum \limits_{j}^{\infty}
\frac{\textstyle 1}{\textstyle  r_{i} r_{j}  }
$
We are facing the problem of an undetermined form
$\infty  \infty$.
Instead of evaluating the potential as a sum over the point charges
we may first rewrite these charges as deltalike charge densities,
$
\rho( r)=
q \sum \limits_{j+}^{\infty} \delta \left( r r_{j+} \right) 
q \sum \limits_{j}^{\infty} \delta \left( r r_{j} \right)
$
and expand these in a Fourier series whose terms determine the Fourier
components $\phi( k)$ of the electrostatic potential.
Since the Fourier representation of a deltafunction requires infinitely
many terms, the Fourier space calculation would again lead to
convergence problems.
The solution is to split up the potential in two wellbehaved
parts, one being represented in $r$space and the other in
$k$space by rapidly converging series.
We demonstrate this on a onedimensional ion lattice with a
charge distribution as depicted in Figure
6.8.
Figure 6.8: Ewald summation
We augment the deltalike point charges by Gaussian charge "clouds"
of opposite sign,
$
\rho'(r) = q_{j}\left(\frac{\textstyle \eta^{2}}{\textstyle \pi} \right)^{3/2}
e^{\textstyle  \eta^{2}( r r_{j})^{2}}
$
(6.11)
to form an auxiliary lattice $1$. A further lattice ($2$) is then introduced
to compensate the additional Gaussian charges, such that
"lattice 1 + lattice 2 = original lattice".
The contributions of the two lattices to the potential are computed
separately:
 Lattice $1$:
Seen from a greater distance, a Gaussian charge cloud resembles a
deltalike point charge, effectively compensating the original charge
it accompanies. The effect of lattice $1$ is therefore best computed in
$r$space, where the series will converge quite rapidly. The
convergence will be faster if the Gaussians are narrow, i.e. if the
parameter $\eta$ in 6.11 is large.
 Lattice $2$:
The potential sum is evaluated in $k$space. When
the Gaussians are broad, i.e. when $\eta$ is small, we will need a smaller
number of Fourier components.
By suitably adjusting $\eta$, optimal convergence of both series may be
achieved.
Proceeding to threedimensional model systems, we consider
a cubic base cell with side length $L$ containing $N$ charges.
Fourier vectors:
$
k \equiv \frac{\textstyle 2 \pi}{\textstyle L} \left( k_{x}, k_{y}, k_{z}\right)
$
with integer $k_{x}$ etc.
Interparticle vectors: including all periodic images of the base
cell, we have
$
r_{i,j, n} \equiv r_{j} + n L  r_{i}
\;\;\; (i,j = 1, \dots, N)
$
where $nL$ is a translation vector in the periodic lattice.
Ewald sum:
$
\phi \left(r_{i}\right) = \frac{\textstyle 4 \pi}{\textstyle L^{3}}
\sum \limits_{j}^{N} q_{j} \left[
\sum \limits_{ k} e^{\textstyle i k \cdot r_{ij}}
k^{2} e^{\textstyle k^{2}/ 4 \eta^{2}} +
\sum \limits_{ n} F(\eta  r_{i,j, n} ) \right]
$
with
$
F(z) \equiv \frac{\textstyle 2}{\textstyle \sqrt{\pi}}
\int_{z}^{\infty} e^{\textstyle t^{2}} dt
$
Note:
Two details need attention:
 The Gaussian charge clouds will formally interact with themselves,
giving rise to a spurious contribution to the potential energy;
this contribution must be subtracted in the final formula.
 The consistent way of taking the infinitesize limit is the following:
 consider a finite (roughly spherical) array of image
cells; surround them by a continuum with some arbitrary dielectric constant
$\epsilon_{s}$, which is usually taken to be $=1;$
 take the limit of an infinitely large repeated array;
this limit still contains a contribution from $\epsilon_{s}$.
Considering these two corrections, we have for the total
potential energy
$
E_{pot} = \frac{\textstyle 1}{\textstyle 2}
\sum \limits_{i}^{N} q_{i} \phi \left( r_{i}\right)
 \frac{\textstyle \eta}{\textstyle \sqrt{\pi}}
\sum \limits_{i}^{N}q_{i}^{2}
+\frac{\textstyle 2 \pi}{\textstyle 3L^{3}}
 \sum \limits_{i}^{N} q_{i} r_{i} ^{2}
$
Another interesting class of particles are those with embedded
point dipoles. Several methods have been devised to deal with
the long range contributions in these model systems. One is a modification
of the Ewald sum; it is known as the
"EwaldKornfeld summation" technique. Other strategies are
the reaction field method and Ladd's
multipole expansion method; see [
VESELY 78] and [
ALLEN 90].
6.5.2 ParticleMesh Methods (PM and P3M):
Huge model systems with $1/r$ potentials (celestial masses or microscopic
charged particles) may be described by introducing
"superparticles" consisting of some $10^{4}10^{8}$
ions, electrons, or stars.[
HOCKNEY 81]
Neglecting the shortrange interactions, Hockney suggested a method
to speed up the dynamics due to the farreaching $1/r$ potential.
Consider a square cell of $M \times M$
subcells of side length $\Delta x = \Delta y = \Delta l$=.
Each subcell should contain an average of $10100$ superparticles.
Equation of motion for superparticle $k$:
$
\ddot{ r}_{k} = \frac{\textstyle q_{k}}{\textstyle m_{k}} \nabla_{k}\Phi =
\frac{\textstyle q_{k}}{\textstyle m_{k}} E( r_{k})
$
where $\Phi( r)$ is the electrostatic or gravitational potential.
It is determined by the charge (or mass) density $\rho( r,t)$
defined by the positions of all superparticles.
To compute $\Phi( r)$ at time $t_{n}$ at the centers of the
subcells, the given configuration of superions is first replaced by a
latticelike charge distribution $\rho_{i,j}$. The easiest
discretization method is the nearest grid point (NGP) rule:
$
\rho_{i,j} = \frac{\textstyle 1}{\textstyle (\Delta l)^{2}}
\sum \limits_{k}^{N} q_{k}
\;\delta \left( \frac{\textstyle x_{k}}{\textstyle \Delta l}i \right)
\delta \left( \frac{\textstyle y_{k}}{\textstyle \Delta l}j \right)
$
A more refined method of charge assignment than the NGP rule is
the cloud in cell (CIC) prescription:
Appropriate fractions of each charge are distributed to the four (2D) or
eight (3D) nearest cell centers.
These fractions, or weights, are assigned in proportion to the overlap
areas of a square of side length $\Delta l$, centered around the particle
under consideration, and the respective neighbor cells
(see Fig. 6.9).
Figure 6.9: Area weighting according to the CIC (cloudincell) rule
The next step is the calculation of the potential produced by the
charge lattice. Any of the relaxation method of Chapter 5 may be
applied, but the FACR technique as developed by Hockney
is preferred.[
HOCKNEY 81]
Result: values of the potential $\Phi_{i,j}$ at the cell
centers. The field strength at the position $r_{k}$ of some
superparticle $k$ in cell $(i,j)$ is then
$
\begin{eqnarray}
E_{x}&=&\left[ \Phi_{i+1,j}  \Phi_{i1,j} \right]\left/ 2\Delta l \right.
\\
E_{y}&=&\left[ \Phi_{i,j+1}  \Phi_{i,j1} \right]\left/ 2\Delta l \right.
\end{eqnarray}
$
Next we integrate the equation of motion:
$
\begin{eqnarray}
r_{k}^{n+1}&=&2 r_{k}^{n} r_{k}^{n1}
+\frac{\textstyle q_{k}}{\textstyle m_{k}} (\Delta t)^{2} E_{i,j}^{n} \\
v_{k}^{n}&=&\left[ r_{k}^{n+1} r_{k}^{n1}\right] \left/
2\Delta t \right.
\end{eqnarray}
$
which completes the time step. Here is the prescription once more:
Figure 6.10: Particlemesh method
Particlemesh method:
At time $t_{n}$ the spatial distribution $\{ r_{k} \}$ of the charged
(or gravitating) superparticles is given.

Assign charge densities $\rho_{i,j}$ to the centers of the cells,
either according to the NGP rule
$
\rho_{i,j} =\frac{\textstyle 1}{\textstyle (\Delta l)^{2}}
\sum \limits_{k}^{N} q_{k}
\;\delta \left( \frac{\textstyle x_{k}}{\textstyle \Delta l}i \right)
\delta \left( \frac{\textstyle y_{k}}{\textstyle \Delta l}j \right)
$
or by some more refined method such as CIC
(see Fig. 6.9.).

Compute the potential at the cell centers, preferably by the FACR
method. For the local field within cell $(i,j)$ use the approximation
$
E_{x}=\left[ \Phi_{i+1,j}  \Phi_{i1,j} \right]\left/ 2\Delta l \right.
$
etc.

Integrate the dynamical equations up to $t_{n+1}$, for instance by
using the Verlet scheme
$
r_{k}^{n+1} = 2 \, r_{k}^{n} r_{k}^{n1}
+\frac{\textstyle q_{k}}{\textstyle m_{k}} (\Delta t)^{2} E_{i,j}^{n}
$

The PM technique considers only the action of the total field
by the distant superparticles. If the shortranged forces may
not be neglected, as in the simulation of molten salts, the
Born, Huggins and Mayer potential is included (see
Table 6.1):
$
U(r) = \frac{\textstyle q_{i}q_{j}}{\textstyle 4\pi \epsilon_{0}r}+B_{ij}
e^{\textstyle \alpha_{ij}r}  \frac{\textstyle C_{ij}}{\textstyle r^{6}}
 \frac{\textstyle D_{ij}}{\textstyle r^{8}}
$
Combining the PM method and the molecular dynamics technique [
HOCKNEY 81],
we may take into account the shortranged forces up to a certain interparticle
distance, while the longranged contributions are included
by the particlemesh procedure. This combination of particleparticle and
particlemesh methods has come to be called PPPM or P^{3}M
technique.
vesely 2006
