Algorithms for Total Energy and Forces in Condensed-Matter 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) + 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)

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 no yes static run energy converged ? 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 orH → H’’ = (L†L)-1(H-E1)C:= L†L ~ 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 → Tr[H ] 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 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

14. 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 )

15. 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’

16. 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 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 • Pulay’s 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 = • DFIBS=0 holds for pure plane-wave basis sets (methods 3,6), not for 1,2,3,5.

20. Forces in LAPW

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

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