Concepts and Math Problems in Electronic Structure Calculations - - PowerPoint PPT Presentation

Concepts and Math Problems in Electronic Structure Calculations

Lin-Wang Wang Scientific Computing Group

• Many-body Schrodinger’s equations
• Density functional theory and single particle equation
• Selfconsistent calculation/nonlinear equation/optimization
• Optical properties
• Basis functions for wavefunctions
• Pseudopotentials
• Technical points in planewave calculations

Many body Schrodinger’s equation Schrodinger’s equation (1930’s): the great result of reductionism !

) , ,.. ( ) , ,.. ( } | | | | 1 2 1 {

1 1 , , 2

t r r t i t r r R r Z r r

N N R i i j i j i i i

Ψ ∂ ∂ = Ψ − + − + ∇ −

∑ ∑ ∑

All the material science and chemistry is included in this equation ! The challenge: to solve this equation for complex real systems.

) ,.. ( ) , ,.. (

1 1 N t i N

r r e t r r Ψ = Ψ

− ω

For stationary solution:

) ,.. ( ) ,.. ( } | | | | 1 2 1 {

1 1 , , 2 N N R i i j i j i i i

r r E r r R r Z r r Ψ = Ψ − + − + ∇ −

∑ ∑ ∑

The famous Einstein formula: E=ħω Ground state: the lowest E state; Excited state: higher E state.

Many body wavefunctions Electrons are elementary particles, two electrons are indistinguishable

) ... ... ,.. ( ) ... ... ,.. (

1 1 N i j N j i

r r r r r r r r Ψ = Ψ α

1

2 =

α 1 = α

, Boson: phonon, photon, W-boson, Higgs-boson, …. (usually particles which transmit forces)

1 − = α

, Fermion: electron, proton, neutron, quark,muon, …. (usually particles which constitute the matter) For our case: electron

) ... ... ,.. ( ) ... ... ,.. (

1 1 N i j N j i

r r r r r r r r Ψ − = Ψ

Many body wavefunctions

) ... ... ,.. ( ) ... ... ,.. (

1 1 N i j N j i

r r r r r r r r Ψ − = Ψ

One example of the antisymmetric wavefunction: Slater determinate

N dr dr r r

N N

= Ψ

∫∫∫

... | ) ,... ( |

1 2 1

normalized antisymmetric Φ1(r1) ΦN(rN) Φ1(rN) ΦN(r1) ……… ……………………… ………

= Ψ ) .... ( 1

N

r r

This is the exact solution for:

) ,.. ( ) ,.. ( } | | | | 1 2 1 {

1 1 , , 2 N N R i i j i j i i i

r r E r r R r Z r r Ψ = Ψ − + − + ∇ −

∑ ∑ ∑

The partial differential equation becomes separable

Another way to look at it: variational methods

N N R i i j i j i i i N

dr dr r r R r Z r r r r E .. ) ,.. ( } | | | | 1 2 1 ){ ,.. (

1 1 , , 2 1

Ψ − + − + ∇ − Ψ =

∑ ∑ ∑ ∫∫∫

The ground state corresponds to the optimized state Ψ which is antisymmetric and normalized. So, we can try variational Ψ for whatever expressions we like variational linear eigen value problem nonlinear problem

• n simplified functions

approximation Plug in the Slater determinate for Ψ, we have (Hartree-Fock equation):

) ( ' ) ' ( ) ' ( | ' | ) ( ) ( } ' | ' | ) ' ( | | 2 1 {

2

r E dr r r r r r r dr r r r R r Z

i i j i j j i R

ϕ ϕ ϕ ϕ ϕ ρ = − + − + − + ∇ −

∑∫ ∫ ∑

Some concepts and terminologies Φ1(r1) ΦN(rN) Φ1(rN) ΦN(r1) ……… ……………………… ……… N+1 1 2 N N+2

= Ψ ) .... ( 1

N

r r

Ei Φj(r) : single particle orbital One orbital can only have one electron (2 include spin)

• ------ Pauli exclusion principle

the N occupied single particle orbitals Φ1 Φ2 ΦN , ,… We also have: the unoccupied orbitals ΦN+1 ΦN+2 , ,… Using one of ΦN+1 ΦN+2 , ,… to replace one of Φ1 Φ2 ΦN , ,… the resulting Slater determinant will correspond to one excited state (band gap)

N N ground excited

E E E E − ≈ −

+1

For the lowest excited state:

Energy breakup

dr dr r r r r dr r R r Z dr r r E

R i i i HF tot

' | ' | ) ( ) ' ( 2 1 ) ( | | ) ( ) ( 2 1

2

∫ ∑ ∫ ∫ ∑

− + − + ∇ − = ρ ρ ρ ϕ ϕ

kinetic Electron-ion Electron Coulomb

dr dr r r r r r r

j i i j i j

' | ' | ) ' ( ) ' ( ) ( ) (

,

∑∫

− + ϕ ϕ ϕ ϕ

Exchange energy

Ecorr=Eexact-EHF

Whatever left from HF Coulomb: ~ 40 eV/atom Kinetic: ~ 40 eV/atom Exchange: ~ 20eV/atom Correlation: ~ 4 eV/atom Typical chemical bond: ~ 2 eV Every term is important For chemical accuracy, we need: ~ 0.05 eV/atom

Different configurations: CI 1 2 N N+2 N+1 hole electron Φ1(r1) ΦN(rN) Φ1(rN) ΦN(r1) … …………………………… … Φj,c(r1) Φj,c(rN) … … j,c SDconf(r1,..rN)= CI: configuration interaction i,v

) ,..., ( ) ( ) ,... (

1 1 N config config N

r r SD config C r r

∑

= Ψ

The number of configuration is exponential, only feasible for a few atom systems. Judicious selection of configurations: MP2, coupled-cluster, etc Traditional quantum chemistry approaches

More on variational many-body wavefunctions e e e e e One electron at r will repulse other electrons near r due to Coulomb inter. Correlation effects: Φ1(r1) ΦN(rN) Φ1(rN) ΦN(r1) ……… ……………………… ………

|)] (| ) ( exp[ ) .... ( 1

∑ ∑

− − − = Ψ

ij j i i i N

r r u r r r χ

Jastrow factor Unfortunately, cannot break down the following integration.

N N R i i j i j i i i N

dr dr r r R r Z r r r r E .. ) ,.. ( } | | | | 1 2 1 ){ ,.. (

1 1 , , 2 1

Ψ − + − + ∇ − Ψ =

∑ ∑ ∑ ∫∫∫

Using Monte-Carlo method to do the integration: variational quantum MC.

Diffusion quantum Monte-Carlo approach

) , ,.. ( ) , ,.. ( } | | | | 1 2 1 {

1 1 , , 2

t r r t i t r r R r Z r r

N N R i i j i j i i i

Ψ ∂ ∂ = Ψ − + − + ∇ −

∑ ∑ ∑

This looks like a classical diffusion equation with finite temperature

{ }

) , ( ) , ( ) (

2

t r S t t r S r V D r r r ∂ ∂ = − + ∇ µ

S(r,t) particle density Using classical Monte-Carlo to simulate the random movements of particles in a 3N dimension space. Problem: S is always positive, but ψ has both positive and negative due to antisymmetry the famous sign problem ! Φ1(r1) ΦN(rN) Φ1(rN) ΦN(r1) ……… ……………………… ……… to divide the 3N space into positive and negative compartments, move articles within. Fix nodal approx: use

Another approach: the density matrix method

N N R i i j i j i i i N

dr dr r r R r Z r r r r E .. ) ,.. ( } | | | | 1 2 1 ){ ,.. (

1 1 , , 2 1

Ψ − + − + ∇ − Ψ =

∑ ∑ ∑ ∫∫∫

' ' ) ' , ; ' , ( } | | | | 1 ){ ' ( ) ' (

2 2 1 1 2 2 1 1 1 2 1 2 1 2 2 1 1

dr dr dr dr r r r r R r Z r r r r r r E

R

ρ δ δ

∑ ∫∫∫

− + − + −∇ − − =

N N N

dr dr r r r r r r r r r r r r ... ) ,... , ' , ' ( ) ,... , , ( ) ' , ; ' , (

3 3 2 1 3 2 1 2 2 1 1

Ψ Ψ = ∫∫ ρ

Great, reduce the N variable function into a 4 variable function !! might not be N-representable ! Problem: ρ ) ' , ; ' , (

2 2 1 1

r r r r

• Many necessary conditions to make ρ N-representable
• The ρ is within some hyperdimension convex cone.
• Linear programming optimization approach
• Recent work: Z. Zhao, et.al, it can be very accurate, but it is still

very expensive (a few atoms).

• No known sufficient condition

The density functional theory

N N N

dr dr r r r r r r r r r ... ) ,... , , ( ) ,... , , ( ) (

2 3 2 1 3 2 1 1

Ψ Ψ = ∫∫ ρ

Any single particle ρ(r) is N-representable. Can we use ρ as one basic variable to determine all other things ?

N N i i j i j i i i N

dr dr r r r V r r r r E .. ) ,.. ( )} ( | | 1 2 1 ){ ,.. (

1 1 , 2 1

Ψ + − + ∇ − Ψ =

∑ ∑ ∑ ∫∫∫

V(r) is one basic variable which determines everything. So V ρ , Now, can ρ V ? (ρ uniquely determine V) We need to prove: we cannot have V1 ρ, and V2 ρ.

Density functional theory (continued) We need to prove: we cannot have V1 ρ, and V2 ρ.

N N i i j i j i i i N

dr dr r r r V r r r r E .. ) ,.. ( )} ( | | 1 2 1 ){ ,.. (

1 1 , 2 1

Ψ + − + ∇ − Ψ =

∑ ∑ ∑ ∫∫∫

Suppose this happens, then V1Ψ1 ρ and V2Ψ2 ρ ) , ( ) , (

2 1 1 1

Ψ < Ψ V E V E

• Since Ψ1 is the variational minimum of V1, so:

dr r r V E E dr r r V E E

Coul K Coul K

) ( ) ( ] [ ] [ ) ( ) ( ] [ ] [

1 2 2 1 1 1

ρ ρ

∫ ∫

+ Ψ + Ψ < + Ψ + Ψ Eq(1) ] [ ] [ ] [ ] [

2 2 1 1

Ψ + Ψ < Ψ + Ψ

Coul K Coul K

E E E E

• Since Ψ2 is the variational minimum of V2, so: E

) , ( ) , (

1 2 2 2

Ψ < Ψ V E V Eq(2) ] [ ] [ ] [ ] [

1 1 2 2

Ψ + Ψ < Ψ + Ψ

Coul K Coul K

E E E E Eq(1),(2) contradict with each other, so we cannot have V1 ρ, and V2 ρ We can also prove, smooth ρ is V-representable (i.e, can find a V ρ) In summary, V is a functional of ρ, thus everything is a functional of ρ

Kohn-Sham equation and LDA Ψ[ρ] exists, so:

N N i i j i j i i i N

dr dr r r r V r r r r E .. ) ,.. ( )} ( | | 1 2 1 ){ ,.. (

1 1 , 2 1

Ψ + − + ∇ − Ψ =

∑ ∑ ∑ ∫∫∫

∫ ∫

+ + − + = dr r r V E drdr r r r r E E

xc kin

) ( ) ( ] [ ' | ' | ) ' ( ) ( 2 1 ] [ ] [ ρ ρ ρ ρ ρ ρ

Great, change the problem to a fluid-dynamics like problem, just one func. ρ(r) Problem: DFT proves that Ekin[ρ], Exc[ρ] exist, but they are unknown. Many approx. for Ekin[ρ]: Thomas-Fermi, Gradient Expan., Wang-Teter. approximate Ekin[ρ] by ∫

dr r r

i i i

) ( ) ( 2 1

2ϕ

ϕ ∇ − ∑

∑

=

i i r

r

2

| ) ( | ) ( ϕ ρ , {φi(r)} are orthonormal.

• L. Sham’s idea:

and

Kohn-Sham equation and LDA (continued) Use local density approximation (LDA) for Exc[ρ]:

∫

= dr r E

xc xc

)) ( ( ] [ ρ ε ρ

Find function εxc(ρ) from simple systems: homogeneous electron gas, where The total energy has been calculated by QMC. The Perdew-Zunger paper. Now, we have the LDA formula:

∫ ∫ ∫ ∑∫

+ + − + ∇ − = dr r r V dr r drdr r r r r dr r r E

xc i i i LDA

) ( ) ( )) ( ( ' | ' | ) ' ( ) ( 2 1 ) ( ) ( 2 1

2

ρ ρ ε ρ ρ ϕ ϕ The ground state solution is a minimum of ELDA for variational {φi(r)} The variational minimum condition: (Kohn-Sham equation)

) ( ) ( )} ( 2 1 {

2

r E r r V

i i i LDA

ϕ ϕ = + ∇ −

) ( )) ( ( ' | ' | ) ' ( ) ( r V r dr r r r r V

xc LDA

+ + − = ∫ ρ µ ρ

and:

Selfconsistent calculations

) ( ) ( )} ( 2 1 {

2

r E r r V

i i i

ψ ψ = + ∇ −

N i i ,.., 1

} {

=

ψ

2

| ) ( | ) ( r r

N i i

∑

= ψ ρ

) (r V

Selfconsistency N electron N wave functions

Planewave expansion of the wavefunction

∑

=

q iqr

e q C r ) ( ) ( ψ

Due space representation

) ( ) ( )} ( 2 1 {

2

r E r r V

i i i

ψ ψ = + ∇ −

diagonal in q space diagonal in real space Fast Fourier Transformation between real space ψ(r) and Fourier space C(q).

A parallel Fast Fourier Transformation code

Time for one FFT (sec)

0.3 0.03 0.003 128x128x128 2 8 8 x 2 8 8 x 2 8 8 576x576x576 EPM calc.

• Specially designed for PW elec.

structure calculation.

• Work load balance
• Memory balance
• Minimum communication

FFT

FFT grids

) ( ) ( ' ) ( ) (

2 2 i i grid i

r r V dr r r V ψ ψ

∑ ∫

−

∆Ω =

Gc1<0.5 Gc2, so: It doesn’t have the usual 1/h2 discretization error, it is exact! We are not doing the usual discretization

dq q V e r V

iqr

) ( ) (

∫

∞

= dq q V e r V

Gc iqr

) ( ) ( '

2

∫

=

Note: : , ) ( ) ( '

i i

r V r V ≠

Pseudopotentials The price: need additional nonlocal potential.

' ) ' ( ) ' ( ) ( ) (

, , ,

dr r r W r W r V

ref R ref R ref R nonloc

ϕ ϕ

∫ ∑

= )

KB form:

A few types of eigen value problems (1) Total energy calculations: need all the occupied states (5% of the all lowest eigenstates). Need them inside a

• uter loop

(2) Nonselfconsistent optical property calculations: need a few states at the interior of the spectrum. One shot calculation. (3) Transport problems: a special eigenstate problem, need eigenstates under special boundary conditions.

Total energy problem

) ( ) ( )} ( 2 1 {

2

r E r r V

i i i

ψ ψ = + ∇ −

N i i ,.., 1

} {

=

ψ

Need: (1) The explicit matrix H is only available for very small systems (used in 70’s). (2) For large systems, Hψ is done using FFT (due space representation), so iterative methods are used. (3) Current methods: CG on Grassman’s manifold:

i i H

Min ψ ψ

ij j i

δ ψ ψ =

Under constraint: (4) Band by band algorithm vs all band algorithm

Total energy problem (continue) (1) Davison’s method vs CG method (2) Residual minimization method / direct inversion in the iterative subspace (RMM-DIIS) Using

1 −

=

l l

HR φ

and use {φl} to get the minimum residual

l l l l

H H R ψ ψ ψ ) ( − =

to generate a Krylov subspace {φl} This by itself is very slow, but subspace diagonalization saves the algorithm Doing each band independently, avoid orthogonalization

Total energy problem (continue) (1) Preconditioning: kinetic energy, diagonal precondition (2) Lanczos method: can be very fast. Long Lanczos iteration (10,000) without explicit orth.

Total energy calculation (continue) Lanczos is faster than CG even without precond. Challenge: (1) how to use precond. (2) how to restart.

Total energy calculation (continued) Wish list for total energy calculation algorithms (1) Iterative method based on Hψ (2) Preconditioning, if possible. (3) Restart from previously converged states. (4) Share Krylov space vectors among eigenstates (Lanczos type methods). (5) Avoid frequent orthogonalization among the eigenstates.

Interior eigenstate problem The challenge: H is not explicitly known, cannot be inverted Have to rely on iterative methods. Typically there is a gap in the spectrum, only interested in gap edge states. E index of states

Folded Spectrum Method and Escan Code

i i i

H ψ ε ψ =

i ref i i ref

H ψ ε ε ψ ε

2 2

) ( ) ( − = −

Other methods for interior eigenstates (1) We can try other methods on (H-Eref)2 , e.g, Lanczos (2) Outer / inner loop methods: inner loop try to approximately invert Hy=x. Does it worth it? How do they compare to direct, one-loop method? PN(E) (3) Jacobi-Davison method. E (4) Challenge: current method works on (H-Eref)2 , the condition number is much worse than H. Can any interior eigenstate method be as easy as working

• n H?

(5) Is interior eigenstate problem intrinsically hard for interative methods.

random N i

H P φ ψ ) ( =

Other ideas

) ( ) ( )} ( 2 1 {

2

r E r r V

i i i

ψ ψ = + ∇ −

Using 3D 7 points finite difference formula for ∇2, H is a sparse matrix in real space grid presentation. The resulting H’ can be factorized directly using ~ 200 Ngrid. Then H’y=x can be solved in a linear scaling. This can be used as a preconditioning, or help to solve the original Hy=x. Some problems: there are nonlocal parts in V(r) , thus it is not really diagonal in real space.

The transport problems

) ( ) ( )} ( 2 1 {

2

r E r r V ψ ψ = + ∇ −

We want to solve: for a given E inside a real space domain,

• utside this domain (or at the boundary), we have

special boundary conditions, e.g:

) ) ( exp( ) ) ( exp( ) ( r E ik r E ik r

• −

+

• =

β ψ

) ) ( exp( ) ( r E ik r

• = α

ψ

) ) ( exp( ) ( r E ik r

• =

ψ

) ( ) ( )} ( 2 1 {

2

r E r r V ψ ψ = + ∇ −

) ) ( exp( r E ik

• −

+ β ) ) ( exp( ) ( r E ik r

• = α

ψ