# Algorithms for Total Energy and Forces in Condensed-Matter DFT codes - PowerPoint PPT Presentation

1 / 26

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.

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.

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

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

## 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 v0(r), the ground state energy

is the stationary point of a uniquely defined functional

Kohn & Sham (1966)

[ –2/2m + v0(r) + Veff[r] (r) ] Yj,k(r) = ej,kYj,k(r)

r(r) = j,k | Yj,k( r) |2

in daily practice:

Veff[r] (r) Veff(r(r)) (LDA)

Veff[r] (r) Veff(r(r), r(r) ) (GGA)

### 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

construct new charge density

for all k determine wavefunctions spanning the occupied space

determine occupancies of states

no

yes

static run

energy converged ?

STOP

yes

relaxation run or molecular dynamics

no

forces converged ?

yes

forces small ?

move ions

STOP

yes

no

### DFT methods for Condensed-Matter Systems

• All-electron methods

• Pseudopotential / plane wave method: only valence electrons (which are involved in chemical bonding) are treated explicitly

1) ‘frozen core’ approximation

projector-augmented wave (PAW) method

2) fixed ‘pseudo-wavefunction’ approximation

### 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-atomcorrectly describe the bonding in different environments ? → transferability

### Basis sets used to represent wavefuntions

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

• All-electron: LMTO (atomic orbitals + spherical Bessel function tails, orthogonalized to neighboring atomic centers)

• PAW: plane waves plus projectors on radial grid at atom centers (additive augmentation)

• All-electron or pseudopotential: Gaussian orbitals

• All-electron or pseudopotential: numerical atom-centered orbitals

• pseudopotentials: plane waves

LCAOs

LCAOs

LCAOs

LCAOs

PWs

### Eigenvalue problem: pre-conditioning

• spectral range of H: [Emin, Emax]

in methods using plane-wave basis functions dominated by kinetic energy;

• reducing the spectral range of H: pre-conditioning H → H’ = (L†)-1(H-E1)L-1 orH → H’’ = (L†L)-1(H-E1)C:= L†L ~ H-E1

• diagonal pre-conditioner (Teter et al.)

### 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 ej,k and eigenvectors Yj,k

• practical up to a Hamiltonian matrix size of ~10,000 basis functions

### Eigenvalue problem: iterative methods

• Residual vector

• Davidson / block Davidson methods(WIEN2k option runlapw -it)

• iterative subspace (Krylov space)

• e.g., spanned by joining the set of occupied states {Yj,k} with pre-conditioned sets of residues {C―1(H-E1) Yj,k}

• lowest eigenvectors obtained by diagonalization in the subspace defines new set{Yj,k}

### Eigenvalue problem: variational approach

• Diagonalization problem can be presented as a minimization problem for a quadratic form (the total energy) (1) (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: → minimization of single residua required

### 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 {Ym}:(pre-conditioned) residuum is added to each Ym

• Preserving the occupied subspace (= orthogonalization of residuum to all occupied states):

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

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

lines of fixed r

### Charge self-consistency

Two possible strategies:

• separate loop in the hierarchy (WIEN2K, VASP, ..)

• combined with iterative diagonalization loop (CASTEP, FHImd, …)

‘charge sloshing’

construct new charge density

and potential

construct new charge density

and potential

iterative diagonalization step

of H for fixed r

{Y(i-1)}→ {Y(i)}

|| Y(i) –Y(i-1) ||<d ?

(H-e1)Y<d ?

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

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

No

No

Yes

Yes

No

STOP

No

STOP

### Achieving charge self-consistency

• Residuum:

• Pratt (single-step) mixing:

• Multi-step mixing schemes:

• Broyden mixing schemes: iterative update of Jacobian Jidea: find approximation to c during runtimeWIEN2K: mixer

• Pulay’s residuum minimization

### Total-Energy derivatives

• first derivatives

• Pressure

• stress

• forces

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 positions1st + 3rd term =

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

### Forces in LAPW

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-1)}→ {Y(i)}

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

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

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

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

Forces converged?

Forces converged?

### Car-Parrinello Molecular Dynamics

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

• equations of motion for the wavefunctions and coordinates

• conserved quantity

• in practical application: coupling to thermostat(s)

### Schemes for damped wavefunction dynamics

• Second-order with dampingnumerical 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)

### 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].