Algorithms for total energy and forces in condensed matter dft codes
1 / 26

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

  • Uploaded on
  • Presentation posted in: General

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.

Download Presentation

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

An Image/Link below is provided (as is) to download presentation

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

Presentation Transcript

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)


  • 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



static run

energy converged ?



relaxation run or molecular dynamics


forces converged ?


forces small ?

move ions




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






Implementations, basis set sizes

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

      • 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

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

Two algorithms for self-consistency









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?

Combining DFT with Molecular Dynamics

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


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


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

  • Login