Algorithms for Total Energy and Forces in CondensedMatter DFT codes

1
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

2
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) Veffr (r) Yj,k(r) ej,k
Yj,k(r)
r(r) ?j,k Yj,k( r) 2
in daily practice Veffr (r) ?Veff(r(r))
(LDA) Veffr (r) ?Veff(r(r), ?r(r)
) (GGA)
3
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

4
initialize charge density
initialize wavefunctions
construct new charge density
for all k determine wavefunctions
spanning the occupied space
determine occupancies of states
energy converged ?
yes static run
no
STOP
yes relaxation run or molecular dynamics
no
forces converged ?
yes
forces small ?
move ions
STOP
yes
no
5
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
6
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

7
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
8
Implementations, basis set sizes
9
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
  or H ? H (LL)-1(H-E1) C LL H-E1
• diagonal pre-conditioner (Teter et al.)

10
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

11
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

12
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 ? TrH from eq. (1)
• metals ? minimization of single residua required

13
Algorithms based on the variational principle for
the total energy

• Single-eigenvector methods residuum
  minimization, e.g. by Pulays 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

14
Conjugate-Gradient Method

• Its not always best to follow straight the
  gradient? modified (conjugate) gradient
• one-dimensional mimi-mization of the total
  energy (parameter f j )

15
Charge self-consistency
Two possible strategies

• separate loop in the hierarchy (WIEN2K, VASP, ..)
• combined with iterative diagonalization loop
  (CASTEP, FHImd, )

charge sloshing
16
Two algorithms for self-consistency
No
No
Yes
Yes
No
STOP
No
STOP
17
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
• Pulays residuum minimization

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

19
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
• DFIBS0 holds for pure plane-wave basis sets
  (methods 3,6), not for 1,2,3,5.

20
Forces in LAPW
21
Combining DFT with Molecular Dynamics

• Born-Oppenheimer MD

• Car-Parrinello MD

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

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

24
Comparison of Algorithms (pure plane-waves)
bulk semi-metal (MnAs), SFHIngx code
25
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.

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