38 2. Gauge invariance
Chapter 3
Numerical implementation
Numerical models are essential for predicting material properties and for investigating dy-namics that cannot be determined analytically. Several approaches exist to approximate the dynamics governed by the many-particle Schr¨odinger equation in order to make solu-tions tractable. These models range from time-dependent generalizasolu-tions of static tight-binding models to dynamics of electrons in a pseudopotential over ab initio approaches starting from density functional theory. Optical responses can then be evaluated from many-body perturbation theory [50, 57, 72], by solving Bethe-Salpeter equation [101, 7]or by applying time-dependent density functional theory [90].
This chapter describes an approach based on solving the time-dependent problem in a basis of Kohn-Sham orbitals that are calculated ab initio. This approach has been imple-mented in the publicly available Python package Ultrafast Light-Matter Interaction Code, ulmic.
Ab initio models allow for a systematic approach to determine properties of a large va-riety of materials. For each material under consideration, the calculation of the ground state and the calculation of the optical properties can be separated into two independent steps. Quantum mechanical simulations are then used to determine the real-time response numerically.
Electron dynamics is modeled in a periodic potential using the independent particle approximation. The wave functions are expanded in periodic Bloch functions, and the electronic system is always taken to initially reside in the ground state, i.e. it consists of an ensemble of valence band electrons, which incoherently populates the valence bands.
For illustrative purposes, analytical potentials are also used for one-dimensional and two-dimensional periodic lattices. While the latter approach does not rely on density functional theory, it can be considered equivalent by considering it as the Kohn-Sham potential, and self-consistently determine the external potential that would give rise this particular Kohn-Sham potential.
The first part of the chapter introduces the density functional theory packages that are used in this work. The desired data is extracted from the various codes, and the format of the data is aligned. Afterwards, the implementation of the ulmic-code is described in detail. Lastly, one-dimensional examples are considered using analytical potentials to
40 3. Numerical implementation
demonstrate the robustness of the formalism.
3.1 Integration with density functional theory pack-ages
In order to develop a robust and reliable tool for determining opto-electric properties of a wide range of materials, it is advantageous to leverage the efforts that have been devoted to develop density functional theory codes over the past decades. As many packages exist, it is desirable to develop a flexible interface, such that data from a number of the popular codes is easily extracted and can be used as input for ulmic.
The wave functions contain all relevant information related to the ground state, so ulmiccan only be independent of the choice of exchange-correlation potential if the ground state data is sufficient. When solving the dynamical equations, this translates into assum-ing that the time-dependent density does not change significantly such that the exchange-correlation term is practically time-independent. While density functional theory only allows for determining the ground state, a set of unoccupied excited states can be deter-mined, and these provide a convenient basis when solving the time-dependent problem. The approach pursued here can be considered an approximation to time-dependent density func-tional theory. While TDDFT is capable of including electron-electron interactions and may include time-dependent exchange-correlation terms, it is unclear to what extent approaches based on TDDFT correctly capture the dynamics. Time-dependent exchange-correlation potentials are especially complicated, and calculations based on time-dependent exchange correlation potentials beyond ALDA are rare. Recent calculations with ALDA show that neglecting the time-dependency of the xc-energy has little impact of high-harmonic spectra.
The external density functional theory codes used in this work are Wien2k, GPAW, QUANTUM ESPRESSO and Abinit. Data from the all-electron code Wien2k [14], and pseudopotential-based codes, such as GPAW [32], Abinit [45] and QUANTUM ESPRESSO [43] have been extracted and tested with ulmic. Energies and matrix elements can be obtained readily from GPAW and from Wien2k by using theopticmodule [2]. Data from Abinit and QUANTUM ESPRESSO can be extracted with the help of the code Yambo [79]. The benefit of interfacing the code with multiple existing DFT codes is to make it more accessible to a larger group of researchers and to verify the reproducibility of the various DFT codes. Reproducibility among different DFT codes is of interest, and at least for band gaps, modern codes tend to produce results that are similar up to sufficiently high accuracy [71].
Workflow
Once the self-consistent problem has been solved and the ground state density has been determined, the wave functions for a large number of excited state are calculated. The wave functions are calculated for an equidistant grid of k points, typically a Monkhorst-Pack grid. For each wave function the energyEnk=hψnk|Hˆ|ψnkiis determined, and the matrix
3.1 Integration with density functional theory packages 41
elements of the momentum operator pnmk =hψnk|pˆˆ|ψmki are calculated for pairs of wave functions at each k point. Overlap integrals Snmk,k+∆kα = hψnk|ψmk+∆kαi are calculated between wave functions at neighboring k points. When calculating the overlap integrals, it is important that the gauges for the wave functions have the same periodicity as the Brillouin zone. I.e. for wave functions at kand k+ ∆kα, where both points lie within the first Brillouin zone the overlap is straightforward to calculate:
Snm,k,k+∆kα =hψnk|e−i∆kα·r|ψmk+∆kαi. (3.1) Along the edge of the Brillouin zone, the wave functions at the points lying outside the first Brillouin zone are calculated from the translation of a points in the interior region:
Snm,k,k+∆kα =hψnk|e−i∆kα·r|ψnk+∆kα−GαeiGα·Ri. (3.2) Practically all maintained DFT codes are interfaced with Wannier90, which is a popular program for constructing Wannier functions [83]. The codes are therefore capable of pro-ducing matrices containing the overlap integrals between neighboring wave functions.
At this stage, all of the necessary information for calculating the dynamics in the in-dependent particle approximation can be extracted. I.e. wave functions do not have to be loaded into ulmic, and they can henceforth be discarded.
For further developments, it is possible to calculate the self-energy in order to correct the band gap predicted by e.g. LDA-based calculations. Several options exist to calculate the self-energy. GPAW allows for calculating G0W0 and GW0 with self-consistent energies.
Abinit is capable of the calculating GW-corrected wave functions self-consistently. For Wien2k, the codes exciting, GAP and GAP2 can be used. For QUANTUM ESPRESSO, the codes Yambo [79] and BerkeleyGW [28] can be used.
42 3. Numerical implementation