Algorithms for Total Energy and Forces in Condensed-Matter DFT codes

IPAM workshop “Density-Functional Theory Calculations for Modeling Materials and Bio-Molecular Properties and Functions”

Oct. 31 – Nov. 5, 2005

P. Kratzer

Fritz-Haber-Institut der MPG D-14195 Berlin-Dahlem, Germany

DFT basics

Kohn & Hohenberg (1965) For ground state properties, knowledge of the electronic density r

(r)

is sufficient. For any given external potential

v

0 (

r

)

,

the ground state energy is the stationary point of a uniquely defined functional

E v

0 [ r ]  F [ r ]  

dr

r (

r

)

v

0 (

r

) Kohn & Sham (1966) [ –  2 /2m +

v

0 (

r

) +

V

eff [ r ] (

r

) ] Y j,

k

(

r

) = e j,

k

Y j,

k

(

r

) in daily practice:

V V

eff eff [ [ r r ] ] (

r

) 

V

eff ( (

r

) 

V

eff ( r (

r

) =  j,

k

r (

r

)) (LDA) r (

r

),  r (

r

) ) (GGA) | Y j,

k

(

r

) | 2

Outline

• flow chart of a typical DFT code • basis sets used to solve the Kohn-Sham equations • algorithms for calculating the KS wavefunctions and KS band energies • algorithms for charge self-consistency • algorithms for forces, structural optimization and molecular dynamics

initialize charge density initialize wavefunctions

for all k

construct new charge density determine wavefunctions spanning the occupied space determine occupancies of states no forces converged ?

no yes energy converged ?

yes

static run

STOP yes

relaxation run or molecular dynamics

move ions no forces small ?

yes STOP

DFT methods for Condensed-Matter Systems

• All-electron methods 1) ‘frozen core’ approximation projector-augmented wave (PAW) method 2) fixed ‘pseudo-wavefunction’ approximation • Pseudopotential / plane wave method: only valence electrons (which are involved in chemical bonding) are treated explicitly

Pseudopotentials and -wavefunctions

• idea: construct ‘pseudo-atom’ which has the valence states as its lowest electronic states • preserves scattering properties and total energy differences • removal of orbital nodes makes plane-wave expansion feasible • limitations: Can the pseudo-atom correctly describe the bonding in different environments ? → transferability

Basis sets used to represent wavefuntions

1. All-electron: atomic orbitals + plane waves in interstitial region (matching condition)

LCAOs LCAOs LCAOs LCAOs PWs

2. All-electron: LMTO (atomic orbitals + spherical Bessel function tails, orthogonalized to neighboring atomic centers) 3. PAW: plane waves plus projectors on radial grid at atom centers (additive augmentation) 4. All-electron or pseudopotential: Gaussian orbitals 5. All-electron or pseudopotential: numerical atom-centered orbitals 6. pseudopotentials: plane waves

Implementations, basis set sizes

5 6 1 2 3 4 implementation (examples) WIEN2K TB-LMTO CP-PAW, VASP, abinit Gaussian, Quickstep, … Dmol 3 ~50 bulk compound ~200 ~50 100..200

50-500 CPMD, abinit, sfhingx, FHImd 100..500

surface, oligo-peptide ~20,000 ~1000 5x10 3 …5x10 5 10 3 …10 4 ~1000 10 4 …10 6

Eigenvalue problem: pre-conditioning

• spectral range of H: [E in methods using plane-wave basis functions dominated by kinetic energy;

E

max min ~ , E max

E cut

 ]  2

k

2

cut

( 2

m e

) • reducing the spectral range of H: pre-conditioning

H → H’ = (L † ) -1( H-E1)L -1 or H → H’’ = (L † L) -1 (H-E1)

C

:= L † L ~ H-E1

• diagonal pre-conditioner (Teter et al.)

C

 16

x

4  8

x

3 8

x

3  12 

x

12 2 

x

2 18 

x

18

x

 27  27 ;

x



T

ˆ /

E cut

Eigenvalue problem: ‘direct’ methods

• suitable for bulk systems or methods with atom centered orbitals only • full diagonalization of the Hamiltonian matrix • Householder tri-diagonalization followed by – QL algorithm or – bracketing of selected eigenvalues by Sturmian sequence → all eigenvalues • practical up to a Hamiltonian matrix size of ~10,000 basis functions e j,k and eigenvectors Y j,k

Eigenvalue problem: iterative methods

• Residual vector

R m

 (

H

 e

m

S

) Y

m

• Davidson / block Davidson methods (WIEN2k option runlapw -it) – iterative subspace (Krylov space) – e.g., spanned by joining the set of occupied states { pre-conditioned sets of residues {C ―1 (H-E1) Y j,k } with subspace defines new set { Y j,k } Y j,k } – lowest eigenvectors obtained by diagonalization in the

Eigenvalue problem: variational approach

• Diagonalization problem can be presented as a minimization problem for a quadratic form (the total energy)

E tot

 

j

.

k

Y

j

,

k

|

T



V ion



V Hartree



V xc

| Y

j

,

k



E ion



ion

(1)

E tot



j

 .

k

e

j

,

k

 

E Hartree

 

E Hartree



E ion



ion

(2) • typically applied in the context of very large basis sets (PP-PW) • molecules / insulators: only occupied subspace is required → Tr[H ] from eq. (1) • metals: 

f i

,

k

e

i

,

k

f i

,

k

→ minimization of single residua required   exp( e

i

,  /

k B T

)  1   1

Algorithms based on the variational principle for the total energy

• Single-eigenvector methods: residuum minimization, e.g. by Pulay’s method • Methods propagating an eigenvector system { Y (pre-conditioned) residuum is added to each Y m m }: – Preserving the occupied subspace

(= orthogonalization of residuum to all occupied states):

• conjugate-gradient minimization • ‘line minimization’ of total energy

Additional diagonalization / unitary rotation in the occupied subspace is needed ( for metals ) !

– Not preserving the occupied subspace: • Williams-Soler algorithm • Damped Joannopoulos algorithm

Conjugate-Gradient Method

• It’s not always best to follow straight the gradient → modified (conjugate) gradient 

m



R

(

m

) |

R

(

m

)

R

(

m

 1 ) |

R

(

m

 1 )

d

(

m

) 

R

(

m

)  

m d

(

m

 1 ) • one-dimensional mimi mization of the total energy (parameter f

j

) Y

j

(

i

 1 )  cos f

j

Y

j

(

i

)  sin f

j d j

(

i

)

Charge self-consistency

Two possible strategies: • separate loop in the hierarchy (WIEN2K, VASP, ..) • combined with iterative diagonalization loop (CASTEP, FHImd, …) ‘charge sloshing’ lines of fixed r

Two algorithms for self-consistency

construct new charge density and potential construct new charge density and potential iterative diagonalization step of H for fixed r No ( H e1)Y < d ?

Yes || r (i) – r (i-1) ||= h ?

No { Y (i-1) }→ { Y (i) } No || Y (i) – Y (i-1) ||< d ?

Yes || r (i) – r (i-1) ||= h ?

STOP No STOP

Achieving charge self-consistency

• Residuum:

R

[ r

in

]  

f j

,

k

• Pratt (single-step) mixing:

j

,

k

Y

j

,

k

2  r

in

r (

i

 1 )

in

 r (

i

)

in

 

R

[ r (

i in

) ] • Multi-step mixing schemes: r (

i in

 1 )  r (

i in

) 

i

  

j R

[ r (

in j

) ]

j



i N

– Broyden mixing schemes: iterative update of Jacobian

J

R

[ r ]  

J

 ( r  r

sc

);

J

 1  c  (

V

Hartree

[ idea: find approximation to c during runtime r ] 

V

XC

[ r ]); WIEN2K: mixer – Pulay’s residuum minimization 

j

 

l l

, 

k A lj

 1

A lk

 1 ;

A lk



R

[ r

in

(

k

) ] |

R

[ r

in

(

l

) ]

Total-Energy derivatives

E SCF

(

n

 d

n

) 

E SCF

(

n

) 

O

( d

n

2 ) • first derivatives – Pressure – stress – forces 

p ij F

    

E



E

 e 

V ij

  

E



R

 Formulas for direct implementation available !

• second derivatives – force constant matrix, phonons Extra computational and/or implementation work needed !

Hellmann-Feynman theorem

• Pulay forces vanish if the calculation has reached self consistency and if basis set orthonormality persists independent of the atomic positions

F

  

dE dR

   

j

,

k

d

Y

j

,

k

|

H

| Y

j

,

k

dR

  Y

j

,

k

d

H

|

dR

 | Y

j

,

k

 Y

j

,

k

|

H

|

d

Y

j

,

k

dR

 1 st + 3 rd term = 

F



IBS



j

 ,

k

e

j

,

k

d dR

 Y

j

,

k

| Y

j

,

k

•  F IBS =0 holds for pure plane-wave basis sets (methods 3,6), not for 1,2,3,5.

Forces in LAPW

F



F HF



F core



F IBS



F T



F FAC F



HF



j



k

, Y

j

,

k

| 

V eff



R

 | Y

j

,

k

 

E ion



ion



R



F



core

  

n core

 (

r

) 

V eff

(

r

)

d

r

F



T



j

 ,

k

r MT

 d Y

j

,

k

T

ˆ Y

j

,

k

dS

 

r MT

 d Y

j

,

k

T

ˆ Y

j

,

k

dS F



FAC

 0

Combining DFT with Molecular Dynamics

• Born-Oppenheimer MD • Car-Parrinello MD move ions move ions construct new charge density and potential construct new charge density and potential { Y (i-1) }→ { Y (i) } || Y (i) – Y (i-1) ||=0 ?

{ Y (i-1) }→ { Y (i) } || Y (i) – Y (i-1) ||=0 ?

|| r (i) – r (i-1) ||=0 ?

Forces converged?

|| r (i) – r (i-1) ||=0 ?

Forces converged?

Car-Parrinello Molecular Dynamics

• treat nuclear and atomic coordinates on the same footing: generalized Lagrangian

L



j

 ,

k



j

,

k

j

,

k



E tot

  

j

,

k

,    • equations of motion for the wavefunctions and coordinates 

j

,

k

 

H

Y

j

,

k

  

j

,

n

Y

n

,

k

M



R

  

F

 • conserved quantity

n

 

j

,

k

j

,

k



E tot

  

j

,

k

j

,

k

• in practical application: coupling to thermostat(s) ,   

Schemes for damped wavefunction dynamics

• Second-order with damping 

j

,

k

  (

H

 e

j

,

k

) Y

j

,

k

  Y

j

,

k

numerical solution: integrate diagonal part (in the occupied subspace) analytically, remainder by finite-time step integration scheme (damped Joannopoulos), orthogonalize after advancing all wavefunctions • Dynamics modified to first order (Williams-Soler) (

i

 1 )

j

,

k

j

,

k

  d (

H

 e

j

,

k

) Y

j

,

k

 exp[  d  (

T

ˆ 

V avg

 e

j

,

k

)] Y

j

(

i

,

k

)  d  (

V

ˆ 

V avg

) Y

j

(

i

,

k

)

Comparison of Algorithms (pure plane-waves)

bulk semi-metal (MnAs), SFHIngx code

Summary

• Algorithms for eigensystem calculations: preferred choice depends on basis set size.

• Eigenvalue problem is coupled to charge-consistency problem, hence algorithms inspired by physics considerations.

• Forces (in general: first derivatives) are most easily calculated in a plane-wave basis; other basis sets require the calculations of Pulay corrections.

Literature

• G.K.H. Madsen et al., Phys. Rev. B 64, 195134 (2001) [WIEN2K].

• W. E. Pickett, Comput. Phys. Rep. 9, 117(1989) [pseudopotential approach].

• G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996) [comparison of algorithms].

• M. Payne et al., Rev. Mod. Phys. 64, 1045 (1992) [iterative minimization].

• R. Yu, D. Singh, and H. Krakauer, Phys. Rev. B 43, 6411 (1991) [forces in LAPW]. • D. Singh, Phys. Rev. B 40, 5428(1989) [Davidson in LAPW].