coulomb interactions
play

Coulomb interactions: P 3 M, MMMxD, ELC and ICC Axel Arnold - PowerPoint PPT Presentation

Coulomb interactions: P 3 M, MMMxD, ELC and ICC Axel Arnold http://www.icp.uni-stuttgart.de Institute for Computational Physics Universit at Stuttgart October 10, 2012 Electrostatics Coulomb energy Pair energy summation


  1. Coulomb interactions: P 3 M, MMMxD, ELC and ICC ∗ ∗ ∗ Axel Arnold http://www.icp.uni-stuttgart.de Institute for Computational Physics Universit¨ at Stuttgart October 10, 2012

  2. Electrostatics Coulomb energy Pair energy summation Bjerrum length e 2 N � ′ U = l B q i q j l B = 4 πǫ 0 ǫ r k B T 2 | r ij | http://www.icp.uni-stuttgart.de i , j = 1 summing up 1 / r electrostatic prefactor ∝ Coulomb pair potential inverse temperature Bjerrum length l B for two unit charges: 1k B T 1l B A. Arnold Coulomb interactions 2/26

  3. Electrostatics Coulomb energy Pair energy summation Potential summation N N � � ′ q i q j U = 1 U = l B q i φ ( r i ) 2 2 | r ij | i = 1 http://www.icp.uni-stuttgart.de i , j = 1 summing up 1 / r potential from solving Coulomb pair potential Poisson’s equation Bjerrum length l B N � ∇ 2 φ ( r ) = − 4 π l B q j δ ( r j − r ) j = 1 equivalent approaches A. Arnold Coulomb interactions 2/26

  4. Electrostatics in periodic boundary conditions Coulomb energy Pair energy summation Potential summation N N ∞ � � � � ′ U = 1 U = l B q i q j q i φ per ( r i ) 2 2 | r ij + m L | i = 1 http://www.icp.uni-stuttgart.de m 2 = S i , j = 1 S = 0 conditionally convergent — solve Poisson’s equation summation order important imposing periodic boundaries numerically difficult U not periodic in coordinates r i U is periodic in coordinates r i these two calculate something different! A. Arnold Coulomb interactions 3/26

  5. Where the difference comes from: the dipole term assume summation in periodic shells surrounded by polarizable material of dielectric constant ǫ ∞ ǫ = ǫ ∞ 3 Pair energy summation 2 1 2 vacuum around: ǫ ∞ = 1 ǫ = 1 ǫ = 1 3 1 0 1 3 Potential summation http://www.icp.uni-stuttgart.de 2 1 2 periodic: ǫ ∞ = ∞ ǫ = ǫ ∞ 3 difference to periodic solution is nonperiodic dipole term � � � 2 2 π U ( d ) = q i r i ( 1 + 2 ǫ ∞ ) L 3 i metallic boundary conditions ǫ ∞ = ∞ always safe never use ǫ ∞ < ∞ for conducting systems A. Arnold Coulomb interactions 4/26

  6. Electrostatics in ESPResSo requires myconfig.h-switch ELECTROSTATICS switching on: inter coulomb < l B > <method > <parameters > methods and their parameters: next 2 hours switching off: http://www.icp.uni-stuttgart.de inter coulomb 0 getting l B , method and parameters: inter coulomb returns e. g. {coulomb 1.0 p3m 7.75 8 5 0.1138 0.0} {coulomb epsilon 80.0 n_interpol 32768 mesh_off 0.5 0.5 0.5} A. Arnold Coulomb interactions 5/26

  7. Assigning charges part 0 pos 0 0 0 q 1 part 1 pos 0.5 0 0 q -1.5 Adding a charged plate constraint plate height < h > sigma < σ > http://www.icp.uni-stuttgart.de plate parallel to xy -plane at z = h , charge density σ requires 2D periodicity (nonperiodic in z ) Adding a charged rod < c x > < c y > lambda < λ > constraint rod center rod parallel to z -axis at ( x , y ) = ( c x , c y ) , line charge density λ requires 1D periodicity (periodic in z ) A. Arnold Coulomb interactions 6/26

  8. The Ewald method statischer Gitterpotentiale , Ann. Phys. 369(3):253, 1921 P. P. Ewald, Die Berechnung optischer und elektro- http://www.icp.uni-stuttgart.de P. P. Ewald, 1888 — 1985 Coulomb potential has 2 problems 1. singular at each particle position 2. very slowly decaying Idea: separate the two problems! one smooth potential — Fourier space one short-ranged potential — real space A. Arnold Coulomb interactions 7/26

  9. Ewald: splitting the potential charge distribution N � � ρ = q i δ ( r − r i − n ) n ∈ L Z 3 i = 1 = + http://www.icp.uni-stuttgart.de replace δ by Gaussians of width α − 1 : � α/ √ π � 3 e − α 2 r 2 ρ Gauss ( r ) = δ ( r ) = ρ Gauss ( r ) + [ δ ( r ) − ρ Gauss ( r )] A. Arnold Coulomb interactions 8/26

  10. The Ewald formula U = U ( r ) + U ( k ) + U ( s ) with � � ′ erfc ( α | r ij + m L | ) U ( r ) = l B q i q j real space correction 2 | r ij + m L | m ∈ Z 3 i , j � U ( k ) = l B 4 π http://www.icp.uni-stuttgart.de k 2 e − k 2 / 4 α 2 | � ρ ( k ) | 2 Gaussians in k -space 2 L 3 k � = 0 � U ( s ) = − α l B q 2 √ π Gaussian self interaction i i forces from differentiation F i = − ∂ U ∂ r i ... coming soon to ESPResSo (on GPU) A. Arnold Coulomb interactions 9/26

  11. Mesh-based Ewald methods M. Deserno and C. Holm, JCP 109:7678, 1998 Computer Simulation Using Particles , 1988 R. W. Hockney and J. W. Eastwood, R. W. Hockney http://www.icp.uni-stuttgart.de J. W. Eastwood replace k -space Fourier sum by discrete FFT discrete FT is exact — constant real space cutoff computational order O ( N log N ) most frequently used methods: P 3 M : optimal method PME SPME A. Arnold Coulomb interactions 10/26

  12. Steps of P 3 M 1. { r i , q i } → ρ ( r ) : interpolate charges onto a grid (window functions: cardinal B-splines) 2. ρ ( r ) → � ρ ( k ) : Fourier transform charge distribution 3. ˆ φ ( k ) = ˆ G ( k )ˆ ρ ( k ) : solve Poisson’s equation by multiplication with optimal influence function ˆ G ( k ) http://www.icp.uni-stuttgart.de (in continuum: product of Green’s function 4 π k 2 and Fourier transform of Gaussians e − k 2 / 4 α 2 ) 4. i k ˆ φ ( k ) → ˆ E ( k ) : obtain field by Fourier space differentiation 4. � E ( k ) → E ( r ) : Fourier transform field back 5. E ( r ) → { r i , F i } : interpolate field at position of charges to obtain forces F i = q i E i A. Arnold Coulomb interactions 11/26

  13. Charge assignment q q 1 1 7 8 2 0.8 3 q q 5 4 6 M ( P ) ( x ) 5 0.6 6 7 0.4 q 3 q 4 0.2 q q 2 1 0 a 0 1 2 3 4 5 6 7 http://www.icp.uni-stuttgart.de x interpolate charges onto h -spaced grid N � ρ M ( r p ) = 1 q i W ( p ) ( r p − r i ) h 3 i = 1 W ( p ) ( r ) cardinal B-splines in P 3 M / SPME A. Arnold Coulomb interactions 12/26

  14. Optimal influence function G opt ( k ) = h 6 i k · � 2 m ∈ Z 3 � ( k + 2 π h m ) � R ( k + 2 π W ( p ) h m ) ˆ �� � 2 2 m ∈ Z 3 � ( k + 2 π | k | 2 W ( p ) h m ) aliasing of continuum force http://www.icp.uni-stuttgart.de R ( k ) = − i k 4 π � k 2 e − k 2 / 4 α 2 with differentiation, Green’s function and transform of Gaussian minimizes the rms force error functional � � � � 2 Q [ F ] := 1 h 3 d 3 r 1 d 3 r F ( r ; r 1 ) − R ( r ) h 3 V A. Arnold Coulomb interactions 13/26

  15. Why to control errors � �� ( F exact − F Ewald ) 2 � � N 1 ∆ F 2 rms force error ∆ F = = N i i = 1 10 r max =1, k max =10 r max =2, k max =10 1 r max =1, k max =20 0.1 ∆ F 0.01 0.001 http://www.icp.uni-stuttgart.de 0.0001 1e-05 0 1 2 3 4 5 α optimal α brings orders of magnitude of accuracy at given required accuracy, find fastest cutoffs compare algorithms at the same accuracy A. Arnold Coulomb interactions 14/26

  16. How to: error estimates 10 total error real space estimate 1 k-space estimate 0.1 ∆ F 0.01 0.001 0.0001 0 1 2 3 4 5 α http://www.icp.uni-stuttgart.de Kolafa and Perram: � q 2 � � 2 i − α 2 r 2 ∆ F real ≈ √ � r max L 3 exp max N Hockney and Eastwood: � � q 2 Q [ ˆ G opt ( k )] i ∆ F Fourier ≈ √ L 3 N A. Arnold Coulomb interactions 15/26

  17. P 3 M in ESPResSo (F. Weik, H. Limbach, AA) tune P 3 M for rms force error τ inter coulomb < l B > p3m tune accuracy < τ > \ [r_cut < r max >] [mesh < n M >] [cao < p >] inter coulomb epsilon ǫ ∞ computes optimal α tunes for optimal speed http://www.icp.uni-stuttgart.de r max real space cutoff (0 to retune) n M = L / h mesh size (0 to retune) p charge assignment spline order p (0 to retune) fixing parameters speeds up tuning, if you know the optimal value second command to set ǫ ∞ (defaults to ∞ (“metallic”)) manually set parameters (dangerous!) < l B > p3m < r max > < n M > < p > < α > inter coulomb A. Arnold Coulomb interactions 16/26

  18. Partially periodic systems http://www.icp.uni-stuttgart.de partially p. b. c. for slablike systems (surfaces, thin films) ... or for cylindrical systems (rods, nanopores) dielectric contrasts at interfaces P 3 M cannot be employed straightforwardly A. Arnold Coulomb interactions 17/26

  19. Another approach: MMM2D far formula e − β √ ( x + k ) 2 +( y + l ) 2 + z 2 � φ β ( r ) = ( x + k ) 2 + ( y + l ) 2 + z 2 � k , l ∈ L Z  � �� = 2 � � � β 2 + p 2 ( y + l ) 2 + z 2  e i px K 0 L p ∈ 2 π l ∈ L Z L Z e − √ β 2 + p 2 + q 2 | z | = 2 π � β 2 + p 2 + q 2 e i px e i qy L 2 � http://www.icp.uni-stuttgart.de p , q ∈ 2 π L Z   e f pq | z | = 2 π  + π e i px e i qy + | z | L 2 β − 1 + O β → 0 ( β ) �  L 2 f pq p 2 + q 2 > 0 screened Coulomb interaction in limit of screening length ∞ other formula for z ≈ 0 optimal computation time O ( N 5 / 3 ) , comparable to Ewald analogously for 1d, but then O ( N 2 ) A. Arnold Coulomb interactions 18/26

  20. Dielectric contrasts z water electrode q ∆ t ∆ b ε t ∆ t q x membrane ε m water l z q ∆ b q ε b http://www.icp.uni-stuttgart.de ∆ b ∆ t q water electrode typical two dimensional systems: thin films, slit pores material boundaries = ⇒ dielectric contrast take into account polarization by image charges can be handled by MMM2D A. Arnold Coulomb interactions 19/26

Download Presentation
Download Policy: The content available on the website is offered to you 'AS IS' for your personal information and use only. It cannot be commercialized, licensed, or distributed on other websites without prior consent from the author. To download a presentation, simply click this link. If you encounter any difficulties during the download process, it's possible that the publisher has removed the file from their server.

Recommend


More recommend