Ammm:Mm ewald

From biowiki
Jump to navigation Jump to search

Truncation of long-range interaction

For condense-phase simualtion, Periodic Boundary Conditon is used to avoid the boundary effect.

It is expensive to compute all pairs of interactions.

Solution: cutoff?

Longer cutoff is (may be) more accurate but more expensive.

Electrostatics

Electrostatic interactions die off very slowly. In crystal the interaction can be rather long-ranged. Cutoff is a bad idea.

We typically use Ewald with dipole correction (tin-foil boundary condition)

Ewald

Detailed description can be found in many books and papers. Here we present the discussion found in book by Frenkel and Smit.[1]

Ewald is a rigorous method to treat long range electrostatic interaction.[2] It is as accurate as infinite cutoff, and the computaitonal cost is rughly similiar to ~ 10 Ang cutoff (when Particle mesh Ewals is used).

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle U_{ele}=1/2\sum _{i=1}^{N}q_{i}\phi \left(r_{i}\right)}

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \phi \left(r_{i}\right)=\sum _{j,n}^{'}{\frac {q_{j}}{\left|r_{ij}+nL\right|}}}

Problem: Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle {\frac {1}{r}}} dies off slowly

Solution: Rewrite the potential in two parts Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle {\frac {1}{r}}={\frac {erfc(r)}{r}}+{\frac {erf(r)}{r}}}


The complementary error function, denoted erfc, is defined in terms of the error function:

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle erfc(r)=1-erf(r)={\frac {2}{\sqrt {\pi }}}\int _{r}^{\infty }e^{-t^{2}}\,dt}

Note: erfc(r)/r dies off much faster, i.e. interaction can be ignored within machine precision beyond a short cutoff.

Physically, this effect can be achieved by adding (two equal and opposite) Gaussian charge distributions to screen the point charges:

  • The Gaussian charge distribution leads to a potential associated with erf(r).
  • The screen potential (point charge plus the negative Gaussian) is associated with erfc(r).

Ewald real space contribution

The screen potential, erfc(r) can now be computed using a short cutoff:

we will show later why the potential of a Gaussian charge distribution is erf via solving Poisson equation.


Ewald reciprocal space contribution

What about the Gaussian charge distribution, the potential, erf(r)/r, diverges in the real-space now --- but dies off quickly in the reciprocal (Fourier) space.

Ewald self correction

Note the equation 12.1.17 include a single summation so there is an interaction between the point charge qi and the (negative) Gaussian charge distribution also located at i (r=0). This contribution can be removed.

First we need to solve the potential due to a (positive) Gaussian.

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \rho _{Gauss}(r)=q_{i}(\alpha /\pi )^{\frac {3}{2}}\exp(-\alpha r^{2})} 

 The Dirac delta function is the limit of Gaussians Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \delta _{a}(x)={\frac {1}{\sigma {\sqrt {\pi }}}}\mathrm {e} ^{-x^{2}/\sigma ^{2}}}
 as Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \sigma }
, the width of Gaussian,Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle \to 0.}
 
or Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \delta_a(r) = \frac{1}{\sigma^{3} {\pi}^{3/2}} \mathrm{e}^{-r^2/\sigma^2}}



Here we proved that the potential due to Gaussian charge distribution is erf. 
When alpha -> infinity, Gaussian charge distribution approaches a point charge and the erf approaches 1. 
Therefore, the potential φ (r) approaches Coulomb's law (1/r). Also notice, to turn off "Ewald", we need alpha=0, 
so we are not adding any (Gaussian} charges to screen the point charges

Ewald energy due to point charges

boundary correction

This contribution to the potential energy corresponds to the k = 0 term that we have neglected thus far. It is permissible to ignore this term if the depolarizing field vanishes. This is the case if our periodic system is embedded in a medium with infinite dielectric constant (a conductor), which is what we have assumed throughout. For simulations of ionic systems, it is essential to use such "conducting" boundary conditions; for polar systems, it is merely advantageous. For a discussion of these subtle points, see [3].

Ewald for multipoles

Detailed derivation is in Smith article in CCP5 newsletter [4]






Particle Mesh Ewald

The is the Ewald coefficient controls how much work is done in the real or reciprocal space. There should be a range of gives the same total energy but with different efficiency. The CPU time required for a fully optimized Ewald summation scales with the number of particles as Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle O(N^{3/2})} . the real-space can use a small cutofff (~9 Ang). The calculation of the Fourier part of the Ewald summation scales as N^2 , which makes the Ewald summation inefficient for large system.

Particle mesh method takes advantage of Fast Fourier Transform to solve the reciprocal contribution. The potential and charge are pre-stored on a grid, which are used later in the interpolation to compute the actual energy and force.[5]

Approximate the potential using interpolation (e.g. order p spline) of precomputed phi on the integer grid.

Early interpolation uses Lagrangian [5] but later B spline, which is smoother (continuous derivative)[6]


The PME scales as N log(N).

Typical use:

  • Real space cutoff 9 Ang for point charge
    • For point multipole, we can afford 6.8 Ang
    • longer cutoff smaller Ewald coefficient (less work in the reciprocal space) will lead to slower speed
  • Reciprocal grid size is <1 Ang (smaller more accurate but slower). X Ang cube box, the no. of grid should be at least 1.2X.
    • For example, 30 A cube box should have 36 grid point.
  • Reciprocal vector K should be big enough (depending on the Ewald coefficient)
    • typical Ewald coefficient is 0.45 when 9 Ang real-space cutoff is used. Larger value needed for shorter real-space cutoff.
  • Spline order: 5 minimum

See http://www.ncbi.nlm.nih.gov/pubmed/15267263 for PME for multipole (the one implemented in TINKER).

Fast Multipole Method

FMM is an alternative O(N) method. It is more efficient than PME for an extreme large system. An algorithm that is of order O(N) is the fast multipole method. The multipole method is based on the idea that a group of particles at a large distance can be considered one big cluster, for which it is not necessary to calculate all particle-particle interactions individually. By clustering the system into bigger and bigger groups, the interactions can be approximated. This approach of Appel [7] leads to an order O(N) algorithm [8]. This algorithm was further refined by Barnes and Hut [9]. In the original algorithm of Appel, the clusters were approximated as a single charge. Greengard and Rokhlin [10] developed an algorithm in which the charge distribution in a cluster is approximated by a multipole expansion. Schmidt and Lee extended this method to systems with periodic boundary conditions [11].

Comparison of Ewald and FMM:

  • Sagui and Darden [12]
  • For conceivable system size, PME is more efficient.

Goddard implementation for million atoms [13]:

vdw interactions

If we truncate the potential at a distance rc, the contribution of the tail of the potential u(r) can be estimated (in three dimensions) using

Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle E_{tail}={\frac {N\rho }{2}}\int _{rc}^{\infty }dru(r)4\pi r^{2}}


Vdw interactions dies off faster, especially the repulsion component Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle 1/r^{12}} . The pairwise interactions are cutoff at certain distance (e.g. 12 Ang). It is possible to correct the ignored long-range interactions for isotropic liquids and additive potentials (RDF approaches unit at the cutoff distance).


Cutoff is less "dangers"; in addition, correction can be made for "homogeneous" system.

Recall the RDF (radial distribution function), or g(r), not only characterizes the microscopic structure but also relates to macroscopic ensemble averages of thermodynamics properties via the following relationship: 


Failed to parse (Conversion error. Server ("https://wikimedia.org/api/rest_") reported: "Cannot get mml. Server problem."): {\displaystyle E=(3/2)NK_{b}T+2\pi N\rho \int _{0}^{\infty }r^{2}v\left(r\right)g\left(r\right)dr}

Beyond the cutoff, the g(r)=1, and the "long tail" correction [14]

  1. ↑ Understanding Molecular Simulation From Algorithms to Applications, Frenkel and Smit, acdemic press 2002
  2. ↑ EE Ewald. Die Berechnung optischer und elektrostatischer Gitterpotentiale. Ann. Phys., 64:253-287, 1921.
  3. ↑ M. Neumann. Dipole moment fluctuation formulas in computer simulations of polar systems. Mol. Phys., 50:841-858, 1983.
  4. ↑ W. Smith CCP5, No 46. Point multipoles in the Ewald Summation revisited 1998
  5. ↑ 5.0 5.1 Darden, York and Pedersen, J. Chem. Phys., Vol. 98, No. 12, 15 June 1993
  6. ↑ Essman, JCP 1995, v103 p8577
  7. ↑ A.W. Appel. An efficient program for many-body simulation. SIAM J. Sci. Stat. Comput., 6:85-103, 1985.
  8. ↑ K. Esselink. The order of Appel's algorithm. Inf. Proc. Lett., 41:141-147, 1992.
  9. ↑ J. Barnes and P. Hut. A hierarchical O(N log N) force-calculation algorithm.Nature, 324:446-449, 1986.
  10. ↑ L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. J. Comp. Phys., 73:325-348, 1987.
  11. ↑ K.E. Schmidt and M.A. Lee. Implementing the fast multipole method in three dimensions. J. Stat. Phys., 63:1223-1235, 1991.
  12. ↑ Annu. Rev. Biophys. Biomol. Struct. 1999. 28:155–79
  13. ↑ J. Chem. Phys., Vol. 97, No. 6, 15 September 1992
  14. ↑ Allen & Tildesley Computer Simulation of liquids
Retrieved from "https://biomolmd.org/mw/index.php?title=Ammm:Mm_ewald&oldid=1566"