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 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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
Recommend
More recommend
