On the adaptive finite element analysis of the Kohn-Sham equations - PowerPoint PPT Presentation

On the adaptive finite element analysis of the Kohn-Sham equations Denis Davydov, Toby Young, Paul Steinmann Denis Davydov, LTM, Erlangen, Germany August 2015

Outline 1. Density Functional Theory and its discretisation with FEM 2. Mesh motion approach 3. A posteriori mesh refinement 4. Numerical examples and discussion 5. Conclusions and future work Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 2 /29

Di fg erent Scales time scale ∂ B n B ∂ B d continuum mechanics deal.II (a finite element D ifferential E quations A nalysis L ibrary) Our goal is to atomistic simulation solve quantum mechanics using h-FEM . quantum physics length scale Bangerth, W., Hartmann, R., & Kanschat, G. (2007). deal.II---A general-purpose object-oriented finite element library. ACM Transactions on Mathematical Software, 33(4), 24–es. doi:10.1145/1268776.1268779 Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 3 /29

Kohn-Sham theorem 1. The ground state property of many-electron system is x ∈ R 3 uniquely determined by an electron density ρ ( x ) 2. There exists an energy functional of the system. In the ground state electron density minimises this functional. Born-Oppenheimer approximation: nuclei are fixed in space, electrons are moving in a static external potential generated by nuclei. O R J R I Electronic state is described by a wavefunction ψ ( x ) Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 4 /29

Density Functional Theory The electron density: Total energy of the system X f α | ψ α ( x ) | 2 ρ ( x ) := E 0 = E kin + E Hartree + E ion + E xc + E zz α partial occupancy orthonormal electronic The kinetic energy of non-interacting electrons: number, such that wavefunctions  � Z � 1 X 2 r 2 f α ψ ∗ E kin := ψ α d x X f α ≡ N e α Ω α α The exchange-correlation energy which represents quantum many-body interactions: Z E xc := ρ ( x ) ε xc ( ρ ( x )) d x Ω The Hartree energy (electrostatic interaction between electrons): E Hartree := 1 ρ ( x ) ρ ( x 0 ) Z Z | x − x 0 | d x 0 d x . 2 Ω Ω Electrostatic interaction between electrons and the external potential of nuclei: " # Z Z Z I X E ion := V ion ( x ) ρ ( x ) d x =: ρ ( x ) d x − | x − R I | Ω Ω I Repulsive nuclei-nuclei electrostatic interaction energy: Z I Z J E zz := 1 X | R I − R J | 2 I , J 6 =I Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 5 /29

Kohn-Sham eigenvalue problem stationary conditions " # � 1 2 r 2 + V e ff ( x ; ρ ) ψ α ( x ) = λ α ψ α ( x ) ρ ( x 0 ) | x − x 0 | d x 0 + δ E xc Z I Z X V e ff ( x ) := V ion ( x ) + V Hartree ( x ) + V xc ( x ) ≡ − | x − R I | + δρ Ω I Explicit Approach : compute Hartree potential more efficiently via solving the Poisson equation need to evaluate for every quadrature V ion �r 2 V Hartree ( x ) = 4 πρ ( x ) point. Hence, bad scaling w.r.t. number of atoms. Implicit Approach : recall that 1/r is Green’s functions of the Laplace operator " # �r 2 [ V Hartree ( x ) + V ion ( x )] = 4 π X solution is not in ! H 1 ρ ( x ) � Z I δ ( x � R I ) I Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 6 /29

Kohn-Sham eigenvalue problem (cont.) stationary conditions " # As a result evaluation � 1 2 r 2 + V e ff ( x ; ρ ) ψ α ( x ) = λ α ψ α ( x ) of scales linearly V ion w.r.t. the number of atoms. V e ff ( x ) := V S V L ⇥ ⇤ ion ( x ) + ion ( x ) + V Hartree ( x ) + V xc ( x ) Gaussian Approach : split Coulomb potential into fully-local short-range and smooth long-range parts erf ( | x − R I | / σ ) X V L ion ( x ) := − Z I | x − R I | I This potential corresponds to the Gaussian charge distribution, i.e. � | x � R I | 2 ✓ ◆ Z I X �r 2 V L ion ( x ) ⌘ � 4 π π 3 / 2 σ 3 exp σ 2 I The remaining short-range part is used directly during assembly of the eigenvalue problem Z I 1 − erf ( | x − R I | / σ ) X X V S | x − R I | − V L ion ( x ) := − ion ( x ) = − Z I | x − R I | I I whereas the Hartree potential and the long-range parts are evaluated by solving the Poisson problem " ◆# � | x � R I | 2 ✓ Z I X V Hartree ( x ) + V L �r 2 ⇥ ⇤ ion ( x ) = 4 π ρ ( x ) � π 3 / 2 σ 3 exp σ 2 I Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 7 /29

Discretisation in space Introduce a finite element basis for the wave-functions and the potential fields ψ α i N ψ X X ϕ i N ϕ ψ α ( x ) = i ( x ) ϕ ( x ) = i ( x ) i i X f α | ψ α ( x ) | 2 ρ ( x ) := α solve iteratively until The generalised eigenvalue problem The generalised Poisson problem convergence i.e. split into a sequence of L ij ϕ j = R H j + R V [ K ij + V ij ] ψ α j = λ α M ij ψ α j linear problems j the potential matrix contribution to the RHS from the electron density Z N ψ i ( x ) V e ff ( x ; ϕ , ρ ) N ψ V ij := j ( x ) d x Z R H Ω 4 π N ϕ j := i ( x ) ρ ( x ) d x Ω the mass (or overlap) matrix possibly non-zero contribution to the RHS from nuclei density - R V Z j N ψ i ( x ) N ψ M ij := j ( x ) d x Ω the kinetic matrix the Laplace matrix K ij := 1 Z Z r N ψ i ( x ) · r N ψ j ( x ) d x r N ϕ i ( x ) · r N ϕ L ij := j ( x ) d x 2 Ω Ω Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 8 /29

Outline 1. Density Functional Theory and its discretisation with FEM 2. Mesh motion approach 3. A posteriori mesh refinement 4. Numerical examples and discussion 5. Conclusions and future work Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 9 /29

Mesh motion approach In order to increase the accuracy of numerical solution and avoid singularities of the Coulomb potential, we need to obtain a mesh such that ions are located at vertexes. - find the closest vertex - prescribe it’s motion to nucleus position O - move other nodes as to… R I I V I Hence we face a problem of determining the u ( V I ) mesh motion field according to the u ( x ) specified point-wise constraints. Approach 1 : Laplace equation Approach 1I : small-strain linear elasticity [ c ( x ) C : r sym u ( x )] · r = 0 on Ω , [ c ( x ) r u ( x )] · r = 0 on Ω , u ( x ) = 0 on ∂ Ω , u ( x ) = 0 on ∂ Ω , such that u ( V I ) = R I � V I . such that u ( V I ) = R I � V I , fourth order isotropic 1 initial position elastic tensor with unit c ( x ) = of the vertex min I | x � R I | . bulk and shear moduli. closest to nuclei I. Jasak, H. & Tukovi ć , Z (2006). Automatic Mesh Motion for the Unstructured Finite Volume Method Transactions of Famena, 30, 1-20. Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 10 /29

Mesh motion (scaled Jacobian element quality) • Naive displacement of the closest vertex to the position of nuclei lead to a unacceptably destroyed elements. • Laplace approach was found to lead to a better element quality near nuclei. Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 11 /29

Outline 1. Density Functional Theory and its discretisation with FEM 2. Mesh motion approach 3. A posteriori mesh refinement 4. Numerical examples and discussion 5. Conclusions and future work Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 12 /29

A posteriori mesh refinement Approach 1 : Residual based error-estimator for each eigenpair: ◆ 2 ✓  � 1 � Z Z 2 r 2 + V e ff � λ α ] 2 d a η 2 e, α := h 2 d x + h e [ [ r ψ α · n ] ψ α e Ω e ∂ Ω e X η 2 η 2 e ≡ e, α α Approach 1I : Kelly error estimator applied to electron density or the electrostatic potential . ϕ ρ Z η 2 e ( • ) := h e [ [ r ( • ) · n ] ] d a ∂ Ω e Marking strategy : marks a minimum subset of a triangulation whose M ⊂ T squared error is more than a given fraction of the squared total error "X # X η 2 η 2 θ ≤ e e e ∈ T e ∈ M Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 13 /29

Outline 1. Density Functional Theory and its discretisation with FEM 2. Mesh motion approach 3. A posteriori mesh refinement 4. Numerical examples and discussion 5. Conclusions and future work Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 14 /29

Treatment of the ionic potential Consider hydrogen atom centered at the origin. V e ff ≡ V ion H V L V ion ( x ) ion ( x ) • The 1/r singularity can not be represented exactly by non-enhanced FE trial spaces. • For the Gaussian charge approach a convergence to the analytical solution is evident. Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 15 /29

Treatment of the ionic potential (cont.) Consider hydrogen and helium atoms that are centered at the origin. H He • Each approach converges to the expected solution with subsequent mesh refinement. • The split of the potential into fully local short-range and smooth long-range part is beneficial from the computational points of view. Otherwise scales badly w.r.t. the number of atoms. Denis Davydov, LTM, Erlangen, Germany College Station (August 2015) 16 /29