Ewald summation

Ewald summation

Ewald summation is a method for computing the interaction energies of periodic systems (e.g. crystals), particularly electrostatic energies. Ewald summation is a special case of the Poisson summation formula, replacing the summation of interaction energies in real space with an equivalent summation in Fourier space. The advantage of this approach is the rapid convergence of the Fourier-space summation compared to its real-space equivalent when the real-space interactions are long-range. Because electrostatic energies consist of both short- and long-range interactions, it is maximally efficient to decompose the interaction potential into a short-range component summed in real space and a long-range component summed in Fourier space.

Derivation

Ewald summation rewrites the interaction potential as the sum of two terms:phi(mathbf{r}) stackrel{mathrm{def{=} phi_{sr}(mathbf{r}) + phi_{lr}(mathbf{r}) where phi_{sr}(mathbf{r}) represents the short-range term that sums quickly in real space and phi_{lr}(mathbf{r}) represents the long-range term that sums quickly in Fourier space. The long-ranged part should be finite for all arguments (most notably r=0) but may have any convenient mathematical form, most typically a Gaussian distribution. The method assumes that the short-range part can be summed easily; hence, the problem becomes the summation of the long-range term. Due to the use of the Fourier sum, the method implicitly assumes that the system under study is infinitely periodic (a sensible assumption for the interiors of crystals). One repeating unit of this hypothetical periodic system is called a "unit cell". One such cell is chosen as the "central cell" for reference and the remaining cells are called "images".

The long-range interaction energy is the sum of interaction energies between the charges of a central unit cell and all the charges of the lattice. Hence, it can be represented as a "double" integral over two charge density fields representing the fields of the unit cell and the crystal lattice:E_{lr} = intint dmathbf{r} dmathbf{r}^{prime} ho_{TOT}(mathbf{r}) ho_{uc}(mathbf{r}^{prime}) phi_{lr}(mathbf{r} - mathbf{r}^{prime})

where the unit-cell charge density field ho_{uc}(mathbf{r}) is a sum over the positions mathbf{r}_{k} of the charges q_{k} in the central unit cell

: ho_{uc}(mathbf{r}) stackrel{mathrm{def{=} sum_{mathrm{charges} k} q_{k} delta(mathbf{r} - mathbf{r}_{k})

and the "total" charge density field ho_{TOT}(mathbf{r}) is the same sum over the unit-cell charges q_{k} and their periodic images

: ho_{TOT}(mathbf{r}) stackrel{mathrm{def{=} sum_{n_{1}, n_{2}, n_{3 sum_{mathrm{charges} k} q_{k} delta(mathbf{r} - mathbf{r}_{k} - n_{1} mathbf{a}_{1} - n_{2} mathbf{a}_{2} - n_{3} mathbf{a}_{3})

Here, delta(mathbf{x}) is the Dirac delta function, mathbf{a}_{1}, mathbf{a}_{2} and mathbf{a}_{3} are the lattice vectors and n_{1}, n_{2} and n_{3} range over all integers. The total field ho_{TOT}(mathbf{r}) can be represented as a convolution of ho_{uc}(mathbf{r}) with a "lattice function" L(mathbf{r})

:L(mathbf{r}) stackrel{mathrm{def{=} sum_{n_{1}, n_{2}, n_{3delta(mathbf{r} - n_{1} mathbf{a}_{1} - n_{2} mathbf{a}_{2} - n_{3} mathbf{a}_{3})

Since this is a convolution, the Fourier transformation of ho_{TOT}(mathbf{r}) is a product

: ilde{ ho}_{TOT}(mathbf{k}) = ilde{L}(mathbf{k}) ilde{ ho}_{uc}(mathbf{k})

where the Fourier transform of the lattice function is another sum over delta functions

: ilde{L}(mathbf{k}) = frac{left(2pi ight)^{3{Omega} sum_{m_{1}, m_{2}, m_{3delta(mathbf{k} - m_{1} mathbf{b}_{1} - m_{2} mathbf{b}_{2} - m_{3} mathbf{b}_{3})

where the reciprocal space vectors are defined mathbf{b}_{1} stackrel{mathrm{def{=} mathbf{a}_{2} imes mathbf{a}_{3} / Omega (and cyclic permutations) where Omega stackrel{mathrm{def{=} mathbf{a}_{1} cdot left( mathbf{a}_{2} imes mathbf{a}_{3} ight) is the volume of the central unit cell (if it is geometrically a parallelepiped, which is often but not necessarily the case). Note that both L(mathbf{r}) and ilde{L}(mathbf{k}) are real, even functions.

For brevity, define an effective single-particle potential

:v(mathbf{r}) stackrel{mathrm{def{=} int dmathbf{r}^{prime} ho_{uc}(mathbf{r}^{prime}) phi_{lr}(mathbf{r} - mathbf{r}^{prime})

Since this is also a convolution, the Fourier transformation of the same equation is a product

: ilde{V}(mathbf{k}) stackrel{mathrm{def{=} ilde{ ho}_{uc}(mathbf{k}) ilde{Phi}(mathbf{k})

where the Fourier transform is defined

: ilde{V}(mathbf{k}) = int dmathbf{r} v(mathbf{r}) e^{-imathbf{k} cdot mathbf{r

The energy can now be written as a "single" field integral

:E_{lr} = int dmathbf{r} ho_{TOT}(mathbf{r}) v(mathbf{r})

Using Parseval's theorem, the energy can also be summed in Fourier space

:E_{lr} = int frac{dmathbf{k{left(2pi ight)^{3 ilde{ ho}_{TOT}^{*}(mathbf{k}) ilde{V}(mathbf{k}) = int frac{dmathbf{k{left(2pi ight)^{3 ilde{L}^{*}(mathbf{k}) left| ilde{ ho}_{uc}(mathbf{k}) ight|^{2} ilde{Phi}(mathbf{k}) = frac{1}{Omega} sum_{m_{1}, m_{2}, m_{3 left| ilde{ ho}_{uc}(mathbf{k}) ight|^{2} ilde{Phi}(mathbf{k})

where mathbf{k} = m_{1} mathbf{b}_{1} + m_{2} mathbf{b}_{2} + m_{3} mathbf{b}_{3}in the final summation.

This is the essential result. Once ilde{ ho}_{uc}(mathbf{k}) is calculated, the summation/integration over mathbf{k} is straightforward and should converge quickly. The most common reason for lack of convergence is a poorly defined unit cell, which must be charge neutral to avoid infinite sums.

Particle mesh Ewald (PME) method

Ewald summation was developed as a method of theoretical physics, long before the advent of computers. However, the Ewald method has enjoyed widespread use since the 1970s in computer simulations of particle systems, especially those interacting via an inverse square force law such as gravity or electrostatics. Applications include simulations of plasmas, galaxies and molecules.

As in normal Ewald summation, a generic interaction potential is separated into two termsphi(mathbf{r}) stackrel{mathrm{def{=} phi_{sr}(mathbf{r}) + phi_{lr}(mathbf{r}) - a short-ranged part phi_{sr}(mathbf{r}) that sums quickly in real space and a long-ranged part phi_{lr}(mathbf{r}) that sums quickly in Fourier space. The basic idea of particle mesh Ewald summation is to replace the direct summation of interaction energies between point particles

:E_{TOT} = sum_{i,j} phi(mathbf{r}_{j} - mathbf{r}_{i}) = E_{sr} + E_{lr}

with two summations, a direct sum E_{sr} of the short-ranged potential in real space

:E_{sr} = sum_{i,j} phi_{sr}(mathbf{r}_{j} - mathbf{r}_{i})

(that's the particle part of particle mesh Ewald) and a summation in Fourier space of the long-rangedpart

:E_{lr} = sum_{mathbf{k ilde{Phi}_{lr}(mathbf{k}) left| ilde{ ho}(mathbf{k}) ight|^{2}

where ilde{Phi}_{lr} and ilde{ ho}(mathbf{k}) represent the Fourier transforms of the potential and the charge density (that's the Ewald part). Since both summations converge quickly in their respective spaces (real and Fourier), they may be truncated with little loss of accuracy and great improvement in required computational time. To evaluate the Fourier transform ilde{ ho}(mathbf{k}) of the charge density field efficiently, one uses the Fast Fourier transform, which requires that the density field be evaluated on a discrete lattice in space (that's the mesh part).

Due to the periodicity assumption implicit in Ewald summation, applications of the PME method to physical systems require the imposition of periodic symmetry. Thus, the method is best suited to systems that can be simulated as infinite in spatial extent. In molecular dynamics simulations this is normally accomplished by deliberately constructing a charge-neutral unit cell that can be infinitely "tiled" to form images; however, to properly account for the effects of this approximation, these images are reincorporated back into the original simulation cell. The overall effect is one of a three-dimensional version of the Asteroids game, in which each dimension "wraps around" on itself. This property of the cell is called a periodic boundary condition. To visualize this most clearly, think of a unit cube; the upper face is effectively in contact with the lower face, the right with the left face, and the front with the back face. As a result the unit cell size must be carefully chosen to be large enough to avoid improper motion correlations between two faces "in contact", but still small enough to be computationally feasible. The definition of the cutoff between short- and long-range interactions can also introduce artifacts.

The restriction of the density field to a mesh makes the PME method more efficient for systems with "smooth" variations in density, or continuous potential functions. Localized systems or those with large fluctuations in density may be treated more efficiently with the fast multipole method of Greengard and Rokhlin.

Dipole term

The electrostatic energy of a polar crystal (i.e., a crystal with a net dipole mathbf{p}_{uc} in the unit cell) is conditionally convergent, i.e., depends on the order of the summation. For example, if the dipole-dipole interactions of a central unit cell with unit cells located on an ever-increasing cube, the energy converges to a different value than if the interaction energies had been summed spherically. Roughly speaking, this conditional convergence arises because (1) the number of interacting dipoles on a shell of radius R grows like R^{2}; (2) the strength of a single dipole-dipole interaction falls like frac{1}{R^{3; and (3) the mathematical summation sum_{n=1}^{infty} frac{1}{n} diverges.

This somewhat surprising result can be reconciled with the finite energy of real crystals because such crystals are not infinite, i.e., have a particular boundary. More specifically, the boundary of a polar crystal has an effective surface charge density on its surface sigma = mathbf{P} cdot mathbf{n} where mathbf{n} is the surface normal vector and mathbf{P} represents the net dipole moment per volume. The interaction energy U of the dipole in a central unit cell with that surface charge density can be written

:U = frac{1}{2V_{uc int frac{left( mathbf{p}_{uc}cdot mathbf{r} ight) left( mathbf{p}_{uc} cdot mathbf{n} ight)dS}{r^{3

where mathbf{p}_{uc} and V_{uc} are the net dipole moment and volume of the unit cell, dS is an infinitesimal area on the crystal surface and mathbf{r}is the vector from the central unit cell to the infinitesimal area. This formula results from integrating the energy dU = -mathbf{p}_{uc} cdot mathbf{dE} where dmathbf{E} represents the infinitesimal electric field generated by an infinitesimal surface charge dq stackrel{mathrm{def{=} sigma dS (Coulomb's law)

:dmathbf{E} stackrel{mathrm{def{=} left( frac{-1}{4piepsilon} ight) frac{dq mathbf{r{r^{3 = left( frac{-1}{4piepsilon} ight) frac{sigma dS mathbf{r} }{r^{3The negative sign derives from the definition of mathbf{r}, which points towards the charge, not away from it.

History

The Ewald summation was developed by Paul Peter Ewald in 1921 (see References below) to determine the electrostatic energy (and, hence, the Madelung constant) of ionic crystals.

caling

Generally different Ewald summation methods give different scalings. Direct calculation gives O(N^2), where N is the number of atoms in the system. The PME method gives O(N,mathrm{log},N), where N is the number of atoms in the system.

ee also

* Paul Peter Ewald
* Madelung constant
* Poisson summation formula

References

* Ewald P. (1921) "Die Berechnung optischer und elektrostatischer Gitterpotentiale", "Ann. Phys." 369, 253-287.

* Darden T, Perera L, Li L and Pedersen L. (1999) "New tricks for modelers from the crystallography toolkit: the particle mesh Ewald algorithm and its use in nucleic acid simulations", "Structure" 7, R55-R60.

* Schlick T. (2002). "Molecular Modeling and Simulation: An Interdisciplinary Guide" Springer-Verlag Interdisciplinary Applied Mathematics, Mathematical Biology, Vol. 21. New York, NY.


Wikimedia Foundation. 2010.

Игры ⚽ Поможем сделать НИР

Look at other dictionaries:

  • Paul Peter Ewald — (January 23, 1888 in Berlin, Germany ndash; August 22, 1985 in Ithaca, New York) was a U.S. (German born) crystallographer and physicist a pioneer of X ray diffraction methods. EducationEwald received his early education in the classics at the… …   Wikipedia

  • Poisson summation formula — The Poisson summation formula is an equation relating the coefficients of the Fourier Series of the periodic extension of a function in terms of the values of the function s continuous Fourier transform. The Poisson Summation Formula was… …   Wikipedia

  • Maillage particulier d'Ewald — Sommation d Ewald La sommation d Ewald (ou parfois somme d Ewald) est une méthode de calcul des énergies d interaction de systèmes périodiques (et particulier des cristaux), et tout particulièrement les énergies électrostatiques. La sommation d… …   Wikipédia en Français

  • Sommation d'Ewald — La sommation d Ewald (ou parfois somme d Ewald) est une méthode de calcul des énergies d interaction de systèmes périodiques (et particulier des cristaux), et tout particulièrement les énergies électrostatiques. La sommation d Ewald est un cas… …   Wikipédia en Français

  • Somme d'Ewald — Sommation d Ewald La sommation d Ewald (ou parfois somme d Ewald) est une méthode de calcul des énergies d interaction de systèmes périodiques (et particulier des cristaux), et tout particulièrement les énergies électrostatiques. La sommation d… …   Wikipédia en Français

  • Wolf summation — The Wolf Summation is a method for computing the electrostatic interactions of systems (eg. crystals). This method is generally more computationally efficient than the Ewald summation. The Energy EquationThe total electrostatic interaction energy …   Wikipedia

  • Molecular dynamics — (MD) is a computer simulation of physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a period of time, giving a view of the motion of the atoms. In the most common version, the trajectories of molecules… …   Wikipedia

  • P3M — For the aircraft, see Martin P3M Particle Particle Particle Mesh (P3M) is a Fourier based Ewald summation method[1][2] to calculate potentials in N body simulations.[3][4][5] The potential could be the electrostatic potential among N point… …   Wikipedia

  • Periodic boundary conditions — In molecular dynamics, periodic boundary conditions (PBC) are a set of boundary conditions used to simulate an effectively infinitely tiled system, usually applied to systems consisting of one or more macromolecules in a bath of explicit solvent …   Wikipedia

  • Multipole expansion — A multipole expansion is a mathematical series representing a function that depends on angles usually the two angles on a sphere. These series are useful because they can often be truncated, meaning that only the first few terms need to be… …   Wikipedia

Share the article and excerpts

Direct link
Do a right-click on the link above
and select “Copy Link”