Processing Math: 11%
To print higher-resolution math symbols, click the
Hi-res Fonts for Printing button on the jsMath control panel.
No jsMath TeX fonts found -- using image fonts instead.
These may be slow and might not print well.
Use the jsMath control panel to get additional information.
jsMath Control PanelHide this Message

jsMath
  Franz J. Vesely > CompPhys Tutorial > Selected Applications > Long Ranged Potentials  
 
 





INDEX
1 FINITE DIFFERENCES
2 LINEAR ALGEBRA
3 STOCHASTICS
4 ORDINARY DE
5 PARTIAL DE
6 SIMULATION
7 QUANTUM M.
8 HYDRODYN.
APPENDICES
< >
Part III: Ch. 6




6.5 Particles and Fields

In many cases the interparticle potential may be safely neglected beyond a few particle diameters. Defining a cutoff distance rco, we may estimate the importance of the neglected "tail" of a potential u(r) by

rcou(r)4r2dr

If the decay of u(r) is slow, this integral may be non-negligible. 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 Particle-Mesh methods are appropriate.

Subsections

6.5.1 Ewald summation

Let
uqq=rq1q2

be the ion-ion interaction between charged particles. In the Ewald summation approach the basic cell containing N2 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 rj+ 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 delta-like 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 delta-function requires infinitely many terms, the Fourier space calculation would again lead to convergence problems.

The solution is to split up the potential in two well-behaved parts, one being represented in r-space and the other in k-space by rapidly converging series. We demonstrate this on a one-dimensional ion lattice with a charge distribution as depicted in Figure 6.8.

Figure 6.8: Ewald summation



We augment the delta-like 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 delta-like 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 three-dimensional 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 infinite-size 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 "Ewald-Kornfeld summation" technique. Other strategies are the reaction field method and Ladd's multipole expansion method; see [ VESELY 78] and [ ALLEN 90].



6.5.2 Particle-Mesh 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 short-range interactions, Hockney suggested a method to speed up the dynamics due to the far-reaching 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 10-100 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 lattice-like 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 (cloud-in-cell) 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_{i-1,j} \right]\left/ 2\Delta l \right. \\ E_{y}&=&-\left[ \Phi_{i,j+1} - \Phi_{i,j-1} \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}^{n-1} +\frac{\textstyle q_{k}}{\textstyle m_{k}} (\Delta t)^{2} E_{i,j}^{n} \\ v_{k}^{n}&=&\left[ r_{k}^{n+1}- r_{k}^{n-1}\right] \left/ 2\Delta t \right. \end{eqnarray}

which completes the time step. Here is the prescription once more:

Figure 6.10: Particle-mesh method

Particle-mesh method:

At time t_{n} the spatial distribution \{ r_{k} \} of the charged (or gravitating) superparticles is given.
  1. 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.).

  2. 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_{i-1,j} \right]\left/ 2\Delta l \right. etc.

  3. 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}^{n-1} +\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 short-ranged 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 short-ranged forces up to a certain interparticle distance, while the long-ranged contributions are included by the particle-mesh procedure. This combination of particle-particle and particle-mesh methods has come to be called PPPM- or P3M technique.


vesely 2006

< >