- 38 Views
- Uploaded on
- Presentation posted in: General

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

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

no

yes

static run

energy converged ?

STOP

yes

relaxation run or molecular dynamics

no

forces converged ?

yes

forces small ?

move ions

STOP

yes

no

- 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

- 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

- 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

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

- 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

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

- 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

- 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

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

- 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

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

- 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

- first derivatives
- Pressure
- stress
- forces
Formulas for direct implementation available !

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

- force constant matrix, phonons

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

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?

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

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

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