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

August 2015 Denis Davydov, LTM, Erlangen, Germany

Denis Davydov, Toby Young, Paul Steinmann

On the adaptive finite element analysis of the Kohn-Sham equations

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/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

2

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Difgerent Scales

Our goal is to solve quantum mechanics using h-FEM. time scale length scale quantum physics atomistic simulation continuum mechanics

∂Bn ∂Bd

B deal.II

(a finite element Differential Equations Analysis Library) Bangerth, W., Hartmann, R., & Kanschat, G. (2007). deal.II---A general-purpose object-oriented finite element library. ACM Transactions

• n Mathematical Software, 33(4), 24–es. doi:10.1145/1268776.1268779

3

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Kohn-Sham theorem

4

• 1. The ground state property of many-electron system is

uniquely determined by an electron density

• 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.

RI O ρ(x) x ∈ R3

Electronic state is described by a wavefunction ψ(x)

RJ

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Density Functional Theory

5

The electron density:

ρ(x) := X

α

fα|ψα(x)|2

Total energy of the system

E0 = Ekin + EHartree + Eion + Exc + Ezz Exc := Z

Ω

ρ(x)εxc(ρ(x)) dx

The exchange-correlation energy which represents quantum many-body interactions:

EHartree := 1 2 Z

Ω

Z

Ω

ρ(x) ρ(x0) |x − x0| dx0 dx .

The Hartree energy (electrostatic interaction between electrons):

Eion := Z

Ω

Vion(x)ρ(x) dx =: Z

Ω

" − X

I

ZI |x − RI| # ρ(x) dx

Electrostatic interaction between electrons and the external potential of nuclei:

Ezz := 1 2 X

I,J6=I

ZIZJ |RI − RJ|

Repulsive nuclei-nuclei electrostatic interaction energy: The kinetic energy of non-interacting electrons:

Ekin := X

α

Z

Ω

fαψ∗

α

 1 2r2

• ψα dx
• rthonormal electronic

wavefunctions partial occupancy number, such that

X

α

fα ≡ Ne

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Kohn-Sham eigenvalue problem

6

stationary conditions

" 1 2r2 + Veff(x; ρ) # ψα(x) = λαψα(x) r2VHartree(x) = 4πρ(x)

Explicit Approach: compute Hartree potential more efficiently via solving the Poisson equation Implicit Approach: recall that 1/r is Green’s functions of the Laplace operator

r2 [VHartree(x) + Vion(x)] = 4π " ρ(x) X

I

ZIδ(x RI) #

solution is not in !

H1

need to evaluate for every quadrature

• point. Hence, bad scaling w.r.t. number of atoms.

Vion Veff(x) := Vion(x) + VHartree(x) + Vxc(x) ≡ − X

I

ZI |x − RI| + Z

Ω

ρ(x0) |x − x0| dx0 + δExc δρ

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Kohn-Sham eigenvalue problem (cont.)

7

stationary conditions

" 1 2r2 + Veff(x; ρ) # ψα(x) = λαψα(x)

Gaussian Approach: split Coulomb potential into fully-local short-range and smooth long-range parts

V L

ion(x) := −

X

I

ZI erf (|x − RI|/σ) |x − RI|

This potential corresponds to the Gaussian charge distribution, i.e.

r2 ⇥ VHartree(x) + V L

ion(x)

⇤ = 4π " ρ(x) X

I

ZI π3/2σ3 exp ✓ |x RI|2 σ2 ◆# Veff(x) := V S

ion(x) +

⇥ V L

ion(x) + VHartree(x)

⇤ + Vxc(x) r2V L

ion(x) ⌘ 4π

X

I

ZI π3/2σ3 exp ✓ |x RI|2 σ2 ◆

The remaining short-range part is used directly during assembly of the eigenvalue problem

V S

ion(x) := −

X

I

ZI |x − RI| − V L

ion(x) = −

X

I

ZI 1 − erf (|x − RI|/σ) |x − RI|

whereas the Hartree potential and the long-range parts are evaluated by solving the Poisson problem As a result evaluation

• f scales linearly

w.r.t. the number of atoms.

Vion

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Discretisation in space

8

ψα(x) = X

i

ψαiN ψ

i (x)

ϕ(x) = X

i

ϕiN ϕ

i (x)

Introduce a finite element basis for the wave-functions and the potential fields

[Kij + Vij] ψαj = λα Mijψαj

The generalised eigenvalue problem

Kij := 1 2 Z

Ω

rN ψ

i (x) · rN ψ j (x) dx

the kinetic matrix

Mij := Z

Ω

N ψ

i (x) N ψ j (x) dx

the mass (or overlap) matrix The generalised Poisson problem

Lij ϕj = RH

j + RV j

Lij := Z

Ω

rN ϕ

i (x) · rN ϕ j (x) dx

the Laplace matrix

RH

j :=

Z

Ω

4πN ϕ

i (x) ρ(x) dx

contribution to the RHS from the electron density possibly non-zero contribution to the RHS from nuclei density - RV

j solve iteratively until convergence i.e. split into a sequence of linear problems

the potential matrix

Vij := Z

Ω

N ψ

i (x)Veff(x; ϕ, ρ)N ψ j (x) dx

ρ(x) := X

α

fα|ψα(x)|2

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/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

9

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Mesh motion approach

10

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.

Jasak, H. & Tuković, Z (2006). Automatic Mesh Motion for the Unstructured Finite Volume Method Transactions of Famena, 30, 1-20.

Hence we face a problem of determining the mesh motion field according to the specified point-wise constraints.

u(x)

Approach 1: Laplace equation

[c(x)ru(x)] · r = 0

• n Ω ,

u(x) = 0

• n ∂Ω ,

such that u(VI) = RI VI , c(x) = 1 minI |x RI| .

initial position

• f the vertex

closest to nuclei I.

Approach 1I: small-strain linear elasticity

[c(x)C : rsymu(x)] · r = 0

• n Ω ,

u(x) = 0

• n ∂Ω ,

such that u(VI) = RI VI .

fourth order isotropic elastic tensor with unit bulk and shear moduli.

• find the closest vertex
• prescribe it’s motion to nucleus

position

• move other nodes as to…

RI VI O I u(VI)

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Mesh motion (scaled Jacobian element quality)

11

• 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)

/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

12

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

A posteriori mesh refinement

13

η2

e,α := h2 e

Z

Ωe

✓  1 2r2 + Veff λα

• ψα

◆2 dx + he Z

∂Ωe

[ [rψα · n] ]2 da

Approach 1: Residual based error-estimator for each eigenpair:

η2

e(•) := he

Z

∂Ωe

[ [r(•) · n] ] da

Approach 1I: Kelly error estimator applied to electron density or the electrostatic potential .

ρ ϕ θ "X

e∈T

η2

e

# ≤ X

e∈M

η2

e

Marking strategy: marks a minimum subset of a triangulation whose squared error is more than a given fraction of the squared total error

M ⊂ T η2

e ≡

X

α

η2

e,α

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/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

14

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Treatment of the ionic potential

15

Consider hydrogen atom centered at the origin.

• 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.

V L

ion(x)

Vion(x) Veff ≡ Vion

H

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Treatment of the ionic potential (cont.)

16

• 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. Consider hydrogen and helium atoms that are centered at the origin.

He H

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Quadrature order

17

For explicit treatment of ionic potential nonlinearities mainly come from the exchange-correlation potential and the Coulomb (1/r) potential. Consider helium and carbon atoms with different quadrature orders used to integrate Hamiltonian matrix and RHS of the Laplace problem.

• For both atoms energy convergence w.r.t. quadrature order is observed

C He

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Quadrature order (cont.)

18

• A posteriori refinement has a bigger influence on the resulting energy in comparison to the

quadrature order used in calculations.

• Insufficient quadrature rule may result in energies lower than those obtained form the radial

solution (i.e. violation of the variational principle).

C He

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Choice of the polynomial degree in FE spaces

19

• Linear FEs are too costly.
• The FE basis in the Poisson problem should be twice the polynomial degree of the eigenvalue

problem in order to maintain variational character of the solution (the difference between energies becomes negative for 2/2 combination).

• The efficiency of the basis to reach the chemical accuracy (1e-3) depends on the problem.

1) which FE basis is the “best” (DoFs per element vs interpolation properties)? 2) can we use the same FE basis for Poisson and eigenvalue problem?

O

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Error estimators

20

Consider neon (Ne) atom.

Ne

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Error estimators (cont.)

21

• Residual error estimator and Kelly estimator applied to eigenvectors result in a good

convergence w.r.t. the expected total energy.

• The error “estimator” based on jumps in density (Kelly) does not lead to convergence.
• The Kelly estimator applied to the solution of the Poisson problem shows suboptimal

convergence. Compare different error estimators:

Ne

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Towards optimisation of nucleus positions

22

Consider a hydrogen molecule ion. Atoms are placed on the x axis symmetrically about the origin with varying interatomic distance. The same structured mesh is used as an input.

H H

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Towards optimisation of nucleus positions

23

Consider carbon monoxide. Atoms are placed on the x axis symmetrically about the origin with varying interatomic distance. The same structured mesh is used as an input.

O C

|ψ7(x)|2

a triple bond

(+) :O ≡ C :(−)

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/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

24

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Summary

25

• Conclusions:

Solution of the quantum many-body problem within the context of DFT using a real- space method based on h-adaptive FE discretisation of the space has been developed.

• Different treatments of the ionic potential produce adequate results given fine

enough meshes and are well behaved with respect to refinement.

• The quadrature rule is not the most influential factor so long as it is accurate

enough to integrate the mass matrix.

• In order to maintain variational nature of the solution, the polynomial degree of the

Poisson FE basis should be twice of that used for the eigenvalue problem.

• The effectiveness of the polynomial degree in the FE basis depends on the problem

considered.

• A residual-based error estimator or Kelly-error estimator applied to eigenvectors

result in adequate convergence. Error estimators based on the density should be avoided.

• The proposed mesh motion approach to move vertices in such a way that each

atom has an associated vertex at its position is well suited to be applied to the structural optimisation of molecules.

Denis Davydov, LTM, Erlangen, Germany College Station (August 2015)

/29

Thank you for your attention

References:

Advanced Grant MOCOPOLY

26 Davydov, D.; Young, T. & Steinmann, P . (2015) On the adaptive finite element analysis of the Kohn-Sham equations: Methods, algorithms, and implementation. Journal for Numerical Methods in Engineering.