Propagation of maximally localized Wannier functions in real-time TDDFT

Yost, Dillon C.; Yao, Yi; Kanai, Yosuke

doi:10.1063/1.5095631

Real-time, time-dependent density functional theory (RT-TDDFT) has gained popularity as a first-principles approach to study a variety of excited-state phenomena such as optical excitations and electronic stopping. Within RT-TDDFT simulations, the gauge freedom of the time-dependent electronic orbitals can be exploited for numerical and scientific convenience while the unitary transformation does not alter physical properties calculated from the quantum dynamics of electrons. Exploiting this gauge freedom, we demonstrate the propagation of maximally localized Wannier functions within RT-TDDFT. We illustrate its great utility through a number of examples including its application to optical excitation in extended systems using the so-called length gauge, interpreting electronic stopping excitation, and simulating electric field-driven quantized charge transport. We implemented the approach within our plane-wave pseudopotential RT-TDDFT module of the QB@LL code, and the performance of the implementation is also discussed.

I. INTRODUCTION

Density functional theory (DFT), based on the Hohenberg-Kohn theorem,¹ is perhaps the most widely used approach to calculate the properties of matter from first principles. However, DFT is limited to ground state properties, and many of the physical, chemical, and biological processes of interest in modern science and technology involve excited state dynamics. Such processes include photoexcitation in dye sensitized solar cells,² hot carrier dynamics in nanomaterials,³ and ionizing radiation in biological materials.⁴ While these phenomena can be studied via advanced spectroscopic methods, first-principles simulations based on time-dependent density functional theory (TDDFT) can be used to provide predictive and detailed insights at the atomistic level.

Currently, time-dependent density functional theory, based on the Runge-Gross theorem,⁵ is one of the most effective and efficient methods for first-principles calculations of excited states and their dynamics.^6,7 Most often, it is the linear response formulation of TDDFT (LR-TDDFT)⁸ that is used, due to its utility in calculating low-energy excitations of molecules and materials, allowing for the prediction and interpretation of electronic excitations and absorption spectra.^9–11 However, as its name suggests, linear response TDDFT can only be applied in the linear response (i.e., weak-field) regime and cannot accurately describe processes involving strong fields, such as laser pulses,¹² and other nonlinear response processes.^13–16 Additionally, even in the linear response regime, LR-TDDFT calculations can come at a prohibitive computational expense if the system includes a large number of electrons¹⁷ or if a broadband spectrum is desired. This is due to the fact that the calculation must be carried out iteratively in a space of occupied and virtual state dimensions.¹⁸ It should be noted, however, that there has been progress in efficient calculations of broadband absorption using the Liouville-Lanczos method¹⁹ and there has been some success in using energy-specific TDDFT to calculate high-energy excited states,²⁰ even though these methods are still limited to the linear response regime.

The real-time propagation approach to time-dependent density functional theory (RT-TDDFT) provides an alternative to LR-TDDFT. Since some of the first uses of real-time propagation approaches in the 1990s, RT-TDDFT has gained popularity for a variety of reasons.²¹ In principle, RT-TDDFT can be used to describe both linear and nonlinear responses of matter to perturbations of any strength. In addition, for large systems and certain properties of interest, the RT-TDDFT approach can be more computationally efficient than LR-TDDFT. A single RT-TDDFT simulation can be used to obtain the entire broadband absorption spectrum,^22,23 even including core excitations.^24–26 Additionally, RT-TDDFT simulations give access to the time-dependent electron density, allowing for molecular-level analysis of the excitation dynamics of interest.^27–29 Due to these factors, in recent years, there has been a surge in applications of RT-TDDFT to a wide range of excited state phenomena such as electronic stopping,^16,30–36 optical absorption,^22,37–39 core electron excitations,^24,26,40,41 electronic circular dichroism spectra,⁴² exciton dynamics in nanostructures,⁴³ atom-cluster collisions,^44,45 and laser-induced water splitting.⁴⁶ The promulgation of RT-TDDFT as a means for simulating excited state phenomena has led to its implementation in a variety of electronic structure codes. These include NWChem,³⁷ Spanish initiative for electronic simulations with thousands of atoms (SIESTA),⁴⁷ CP2K,⁴⁸ scalable ab initio light-matter simulator for optics and nanoscience (SALMON),⁴⁹ Octopus,^50,51 Q-Chem,⁵² GAUSSIAN,^53–55 MOLGW,^34,56 Quantum Espresso,⁵⁷ and QBOX/QB@LL.^58–60 Among these codes, the implementations vary in the underlying basis sets used, with implementations ranging from real-space grids,^47,49,50 Gaussian-type atomic orbitals,^{34,37,47,53,54,56} plane-waves,⁵⁸ and mixed Gaussians/plane waves (PWs).⁴⁸

Recently, several RT-TDDFT implementations have exploited the gauge freedom inherent in the time-dependent Kohn Sham (TDKS) equations that underlie RT-TDDFT calculations. Under arbitrary unitary gauge transformations, the physical properties of the electron dynamics are unchanged, but certain gauge transformations allow for reformulations of the TDKS equations that can be advantageous in certain physical and numerical situations.⁶¹ For instance, in recent studies, Pemmaraju et al.²⁶ and Noda et al.⁴⁹ have implemented velocity-gauge formulations of RT-TDDFT that provide a convenient framework to simulate the responses of periodic systems to both weak and intense electric fields which, in the velocity gauge, can be represented via vector potentials.⁶¹

Another interesting example of gauge transformations in RT-TDDFT, demonstrated in recent studies by Lin and co-workers, is the propagation of electron dynamics in the parallel transport gauge.^62,63 In the parallel transport gauge, the unitary transformation is performed such that oscillations of the propagating orbitals are minimized. This allows for significantly larger propagation time steps to be used with implicit time integration schemes such as the Crank-Nicolson scheme, providing computational advantages in practical RT-TDDFT simulations, especially when hybrid exchange-correlation (XC) functionals are involved.⁶⁴

In this work, we exploit gauge freedom in two regards: we represent the KS orbitals in the gauge of maximally localized Wannier functions (MLWFs) and, as it naturally follows, we represent electric fields in the length gauge corresponding to the formulation of finite fields as scalar potentials. The foundational implementation of RT-TDDFT in the QB@LL branch of the Qbox code^58,59 used in this study involves a plane wave pseudopotential (PP) formalism with periodic boundary conditions (PBC) in which the TDKS states are represented as Bloch orbitals. However, in this work, through computation and application of a unitary transformation matrix at each RT-TDDFT simulation step, we propagate the TDKS orbitals as time-dependent MLWFs (TD-MLWFs). Due to their relationship with the dynamic Berry phase,⁶⁵ the TD-MLWFs allow us to calculate the dynamic polarization throughout RT-TDDFT simulations involving isolated and periodic systems in static or time-dependent electric fields. Through such simulations, we can calculate linear response quantities such as optical absorption spectra and we can also simulate nonlinear responses to strong fields associated with laser pulses or ionizing particle radiation. Additionally, MLWFs are often used to give chemically intuitive representations of the electronic structure of molecules and materials in terms of bonds, lone-pairs, etc., and TD-MLWFs have the same utility in the context of RT-TDDFT. By decomposing the dynamic polarization response in terms of contributions from TD-MLWFs associated with different chemical moieties, we determine, for instance, how a water molecule’s lone-pair electrons contribute to individual peaks in the total optical absorption spectrum. We also use the real-time propagation of MLWFs to study the electron excitation dynamics of a system in response to proton irradiation. Normally, excitations induced by photonic and ionic irradiation are difficult to compare, but our implementation of TD-MLWFs allows us to calculate a type of electronic stopping response spectrum, which can be used to make such comparisons. Additionally, we perform simulations of optically driven quantized charge transport in polyacetylene to demonstrate the length-gauge time-dependent field representation and to consider the possibilities for using TD-MLWF propagation to study dynamic topological phenomena.

II. THEORETICAL AND COMPUTATIONAL DETAILS

A. Time-dependent Kohn Sham equations

The Runge-Gross theorem⁵ allows us to describe the quantum dynamics of a system evolving from a given initial state through a time-evolving electron probability density. This extension of the Hohenberg-Kohn theorem,¹ formulated in the time domain, gave rise to the widely successful time-dependent density functional theory (TDDFT), now frequently used for excited state calculations and simulations. As in ground state DFT, practical TDDFT calculations rely on equations expressed in terms of a system of noninteracting Kohn-Sham (KS) particles in an effective KS potential. The time-dependent Kohn-Sham equations (TDKS) are as follows:

\begin{matrix} i ℏ \frac{d}{d t} | ϕ_{i} (r, t) (⟩ & = & \hat{H} (t) | ϕ_{i} (r, t) (⟩ \\ = & \{\hat{T} + {\hat{V}}_{e x t} (t) + {\hat{V}}_{H X C} [ρ]\} | ϕ_{i} (r, t) (⟩, \end{matrix}

(1)

ρ (r, t) = \sum_{i} {|ϕ_{i} (r, t)|}^{2} .

(2)

In Eq. (1), $\hat{T}$ is the kinetic energy operator $\frac{- ℏ}{2 m} \nabla_{r}^{2}$ and $V_{H X C} [ρ] (r, t) = \int \frac{ρ (r^{'}, t)}{|r - r^{'}|} d r^{'} + \frac{δ E_{X C}}{δ ρ (r, t)}$ is the sum of the Hartree (H) potential and the exchange-correlation (XC) potential, which is derived from a universal XC functional E_XC[ρ]. The electron density ρ here is expressed as the sum of square amplitudes of the occupied KS orbitals ϕ_i (labeled by the state index i). The total energy at time t can be expressed as

\begin{matrix} E (t) & = & \sum_{i} ⟨ϕ_{i} (r, t) |\hat{T}| ϕ_{i} (r, t)⟩ \\ + \int ρ (r, t) V_{e x t} (r, t) d r + E_{H X C} [ρ] (t) . \end{matrix}

(3)

This energy functional of the time-dependent density can be shown to obey energy conservation when the so-called adiabatic approximation⁶⁶ is adopted, which provides a useful observable for validating numerical implementations and practical simulations.⁵⁸ In most RT-TDDFT simulations of realistic systems, the exchange-correlation (XC) potential is approximated with the adiabatic approximation. That is, the XC potential depends only on the instantaneous electron density, neglecting any memory effects. While time-dependent XC functionals beyond the adiabatic approximation are an area of active research,^67–71 such a topic is beyond the scope of this work. Even so, within the adiabatic approximation to the XC functional, TDDFT calculations have still been successful in predicting excited state properties for a diversity of molecular and material systems.

While the majority of applications of TDDFT remain within the linear-response formulation of TDDFT (LR-TDDFT), this real-time propagation approach (RT-TDDFT) has become increasingly popular in recent years. In a RT-TDDFT calculation, the computational workload entails solving a set of time-dependent KS equations, which are nonlinear partial differential equations, in time, for given initial conditions. The time dependence of the external potential in the TDKS system can arise explicitly and implicitly from a variety of sources, depending on the phenomenon of interest. For example, in applications of RT-TDDFT to study electronic stopping processes,^{16,27,30,35,36,72,73} a fast-moving charged ion gives rise to an external potential that varies as the ion’s position changes over the course of time. In addition to ionic motion, the external potential can vary in time in situations involving a dynamic applied electric field. A key advantage of RT-TDDFT is that it provides a framework within which a wide variety of dynamical electronic processes, such as electronic stopping and photoexcitation, can be treated on an equal footing.

B. Bloch and Wannier orbital representations

In the KS and TDKS equations, a Bloch representation form naturally arises in extended systems, in which periodic boundary conditions are adapted.⁷⁴ The TDKS orbitals $φ_{n k} (r, t)$ are written in the Bloch function form

φ_{n, k} (r, t) = u_{n, k} (r, t) e^{i k \cdot r},

(4)

where u_nk(r, t) are Bloch orbitals having the periodicity of the lattice, n is the eigenstate index, and k is the Brillouin zone vector. By preserving the electron probability density from the KS single-particle orbitals, any unitary transformation of the above Bloch functions gives an equivalent physical description of the system. Thus, one can apply a unitary transformation to the Bloch functions that results in a representation in a set of “Wannier functions” (WFs)^75–77 as follows:

w_{n, R} (r, t) = \frac{Ω}{{(2 π)}^{3}} \int_{B Z} d k e^{- i k \cdot R} φ_{n, k} (r, t),

(5)

where the Wannier functions $w_{n R}$ are indexed with a bandlike index n located in the lattice-periodic cell R with a real-space cell volume Ω. There exists a gauge freedom in this procedure due to the fact that one can apply an arbitrary phase transformation $e^{i θ_{n} (k)}$ to the Bloch orbitals without changing the physical observables of the system. However, in terms of the resultant WFs, such a transformation to the Bloch orbitals could alter the shapes and spreads of the corresponding WFs. Thus, due to this gauge freedom, in this most general definition, the WFs are not unique.⁷⁸

In order to address this ambiguity, a variety of schemes have been proposed, including projection methods,^79–81 variational procedures,⁸² and methods based on symmetry considerations.^83,84 However, the most widely used method is the maximally localized Wannier function (MLWF) procedure developed by Marzari and Vanderbilt.⁸⁵ In this method, a spatial spread functional is defined, and then the unitary transformation can be optimized in order to minimize the spread functional, thus maximizing the localization of the WFs. Such maximally localized Wannier functions (MLWFs) typically have the character of being localized on chemical bonds and other familiar chemical moieties, which is why they have gained popularity for their use in interpreting the electronic structures of various matter.^27,78 Throughout the entirety of this work, MLWFs are used. The algorithmic and computational details of the localization procedure are discussed later in the manuscript.

Because the transformation between Bloch and MLWF states is unitary, both representations comprise equally valid descriptions of the electronic band subspace.⁷⁸ Thus, in the TDKS equations, one can substitute the single-particle KS orbitals, normally taken to be the Bloch functions, with MLWFs. The charge density obtained by summing the square amplitudes of the KS orbitals is identical to that obtained by summing the square amplitudes of the MLWFs, meaning that we have an equivalent framing of RT-TDDFT calculations in terms of time-dependent MLWFs (TD-MLWFs). The TDKS equation in the Bloch representation is

\begin{matrix} i \frac{d}{d t} | u_{j, k} (r, t) (⟩ & = & \{) \frac{1}{2 m} {[- i ℏ \nabla_{r} + k]}^{2} + {\hat{V}}_{i o n} (t) \\ + {\hat{V}}_{H X C} [{u_{j, k} (r, t)}] \}) | u_{j, k} (r, t) (⟩, \end{matrix}

(6)

where Eq. (6) can also be written in terms of Wannier functions as

\begin{matrix} i \frac{d}{d t} | w_{l} (r, t) (⟩ & = & \{) {\hat{U}}^{M L} + \frac{1}{2 m} {[- i ℏ \nabla_{r}]}^{2} + {\hat{V}}_{i o n} (t) \\ + {\hat{V}}_{H X C} [{w_{l} (r, t)}] \}) | w_{l} (r, t) (⟩, \end{matrix}

(7)

where $w_{l}$ (r, t) is the time-dependent MLWF (TD-MLWF) with a bandlike index l, where ${\hat{U}}^{M L}$ is the unitary operator that ensures the maximal localization,

{\hat{U}}^{M L} | w_{l} (t) (⟩ \equiv \sum_{m}^{N} U_{l m} (t) | w_{m} (t) (⟩,

(8)

where the unitary transformation matrix U_lm minimizes the spread functional

M i n {\{) \sum_{n}^{N} [) ⟨w_{n} |{\hat{r}}^{2}| w_{n}⟩ - {⟨w_{n} |\hat{r}| w_{n}⟩}^{2}]) \})}_{U},

(9)

where the quantum mechanical (QM) position operator $\hat{r}$ is defined generally even for extended systems⁸⁶

⟨\hat{r}⟩ = \frac{L}{2 π} I m \ln ⟨ψ |e^{\frac{i 2 π}{L} \hat{r}}| ψ⟩

(10)

where $L = ‖ R ‖$ is the cell dimension. The unitary operator in Eq. (7) does not change the quantum dynamics governed by the time-dependent electron density, and therefore it does not affect physical observables. The algorithmic details of the unitary transformation matrix calculation are discussed in Sec. II E.

One of the central motivations for carrying out this transformation to MLWFs is that it allows for the calculation of electric polarization in extended periodic systems. While simple formulas exist for the calculation of electric dipoles of molecular systems with localized wavefunctions, these formulas cannot generally be extended to periodic systems in which the Bloch wavefunctions are delocalized in real space. In extended systems, the analog to the electric dipole moment is the electric polarization, defined as the electric dipole moment per unit volume. While the calculation of this quantity for periodic systems may seem intuitive, the arbitrary choice between different valid unit cells can result in contradictory results for the electric polarization. This infamous problem was elegantly solved by Resta⁸⁷ and King-Smith and Vanderbilt⁸⁸ in the work that is now collectively known as the “modern theory of polarization.”⁸⁹ In the modern theory of polarization, the electric polarization of a periodic system can be formulated in terms of Berry phases, or equivalently, in terms of Wannier functions.

MLWFs have been used in the context of ground-state DFT to calculate the adiabatic stationary resonant state of an insulating material in static electric fields.^90,91 Souza et al. also showed that this concept could be extended to dynamics via the definition of a nonadiabatic Berry phase polarization.⁶⁵ Thus, through the incorporation of MLWFs into RT-TDDFT, we can track the dynamics of a system’s polarization in time. The access to the dynamic polarization of a system in response to an external perturbation allows for the calculation of important quantities such as the transient current and frequency-dependent optical absorption.⁶⁵

C. Finite electric fields and periodic boundary conditions

While it is relatively straight-forward to perform RT-TDDFT calculations with time-varying external potentials caused by ionic motion, as has been done in RT-TDDFT studies of electronic stopping,^30,35 the treatment of spatially homogeneous electric fields requires more careful consideration in cases involving periodic systems. Fundamentally speaking, application of the Runge Gross theorem in cases of an extended periodic system in a homogeneous electric field is not formally justified. Instead, time-dependent current-density functional theory (TDCDFT) should be used due to its incorporation of the macroscopic current.⁹² Despite this formal limitation, RT-TDDFT studies have shown that many response properties can in practice be successfully acquired without TDCDFT. That being said, there are several additional theoretical complications to be considered. One way to treat the response of electrons to electric field excitations in RT-TDDFT is to introduce the time-dependent spatially homogeneous electric field in the KS Hamiltonian through a scalar potential in the so-called “length gauge.” The electric field can be included through an additional potential in the Hamiltonian

{\hat{V}}_{E} (t) = e E (t) \cdot \hat{r},

(11)

where E(t) and $\hat{r}$ are the time-dependent electric field and the quantum-mechanical position operator, respectively. However, the added potential makes the Hamiltonian nonperiodic, as usually required for modeling extended systems with periodic boundary conditions (PBCs). Thus, instead of this length gauge formulation, it is common to move to a different gauge in electromagnetism. The electric field can be equivalently represented as the magnetic flux by the vector potential

A (t) = - c \int^{t} E (t^{'}) d t^{'} .

(12)

In the literature, RT-TDDFT simulations of periodic systems with homogeneous electric field employ the so-called velocity-gauge formulation.^26,41,61,93 The TDKS equation can be written as

\begin{matrix} i ℏ \frac{d}{d t} | u_{j, k} (r, t) (⟩ & = & \{) \frac{1}{2 m} {[- i ℏ \nabla_{r} + k + \frac{e}{c} A (t)]}^{2} \\ + {\hat{V}}_{e x t} (t) + {\hat{V}}_{H X C} [ρ] \}) | u_{j, k} (r, t) (⟩, \end{matrix}

(13)

In the above velocity gauge, the Hamiltonian associated with Eq. (13) is periodic for spatially homogeneous electric fields, allowing for a Bloch wavefunction to be used in the TDKS equations.²⁶ Although the velocity gauge is commonly used when simulating extended periodic systems, it is not without its limitations. In particular, unlike the scalar potential, the vector-potential results in a Hamiltonian which inherently changes nonadiabatically in time, even for the static field case.⁶⁵ Thus, unlike the length gauge, the velocity-gauge is not suitable for calculating the stationary resonant polarization state of an insulator in a static electric field using time-independent DFT.

The use of MLWFs instead of Bloch orbitals allows us to employ the length gauge in which the homogeneous electric field is represented by a scalar potential in the Hamiltonian. For nonmetallic systems with a finite bandgap [termed “Wannier-representable (WR)” systems⁶⁵], the MLWFs are spatially well-localized within each periodic simulation cell and the scalar potential can be applied to each MLWF individually such that the same homogeneous electric field is described. This formalism has already been demonstrated for the static case⁹⁴ in the context of first-principles molecular dynamics (MD) simulations.⁹⁵ Due to their connection to the dynamic Berry phase polarization (discussed in more detail in Sec. II D), MLWFs provide a computationally attractive, and physically intuitive, avenue for simulating systems in finite electric fields with RT-TDDFT in the MLWF/length-gauge. This TDKS/TD-MLWF equation can be written as

\begin{matrix} i ℏ \frac{d}{d t} | w_{l} (r, t) (⟩ & = & \{) {\hat{U}}^{M L} + [) \frac{1}{2 m} {[- i ℏ \nabla_{r}]}^{2} + {\hat{V}}_{e x t} (t) \\ + {\hat{V}}_{H X C} [ρ] + e E (t) \cdot \hat{r}]) \}) | w_{l} (r, t) (⟩ . \end{matrix}

(14)

Including the unitary operator ${\hat{U}}^{M L}$ here, as in Eq. (7), ensures that the Wannier functions remain maximally localized during the propagation, and this allows us to use the length gauge for applying a homogeneous electric field even when periodic boundary conditions are adapted. Examples of RT-TDDFT propagation of TD-MLWFs perturbed by an impulsive electric field are shown in Fig. 1. In practice, preserving the maximal locality of the Wannier functions enables us to easily calculate the positions of the MLWF centers at each time step via the diagonal elements of the quantum-mechanical position operator matrices [discussed in more detail in Sec. II E 1, see Eq. (23)].

FIG. 1.

View large Download slide

Isosurfaces of a single TD-MLWF orbital in an isolated water molecule (upper) and in crystalline silicon (lower). Each panel shows the TD-MLWFs at a different point in time in RT-TDDFT simulations after the systems have been perturbed by an instantaneous electric field impulse.

D. Dynamic polarization via Wannier functions

The properties of insulating crystals in static electric fields can be calculated via the iterative determination of field-dependent WFs through the minimization of the so-called “electric enthalpy” functional.⁹⁰ This formalism was also generalized and extended to the time-dependent domain in the work by Souza et al.,⁶⁵ in which it was shown that the dynamic polarization can be expressed as a nonadiabatic geometric phase. The dynamic current is defined as the rate of polarization change per unit volume with respect to time

J (t) = \frac{d P_{e l} (t)}{d t},

(15)

where the dynamic electronic polarization P_el(t) is given by the valence-band Berry phase^88,90,96

P_{e l} (t) = \frac{- 2 e}{{(2 π)}^{3}} \sum_{n} \int_{B Z} d k ⟨u_{n, k} (r, t) |i \nabla_{k}| u_{n, k} (r, t)⟩,

(16)

where u_n,k are the Bloch functions characterized by a band index n. In this time-dependent case, we do not assume that changes in the Hamiltonian are adiabatic. Instead, Eq. (16) can be interpreted as the nonadiabatic geometric phase.⁶⁵ Alternatively, this Berry phase expression can be transformed into a real-space representation in terms of the occupied Wannier functions such as MLWFs

P_{e l} (t) = - \frac{2 e}{Ω} \sum_{i} ⟨w_{i} (r - R, t) |\hat{r}| w_{i} (r - R, t)⟩,

(17)

where the dynamic polarization is recast in real space as the vector sum of the charge centers of mass of the MLWFs (MLWFCs), $w_{i}$ ⁠, which correspond to the expectation values of the quantum mechanical position operator $\hat{r}$ for a periodic system. Here, the location of the MLWFC is the indeterminate modulus of the lattice vector R, and consequently the dynamic Berry-phase polarization P_el(t) is also indeterminate mod $- \frac{2 e R}{Ω}$ ⁠. This uncertainty is the so-called “quantum of polarization.” Thus, as described in the modern theory of polarization, it is actually the change in polarization that is the physical quantity that needs to be measured.

Using MLWFs, one can obtain this physically intuitive definition of the polarization, which is expressed in terms of the geometric centers of charge of the MLWFs (often called “Wannier centers,” or WCs). As illustrated by Eq. (17), the dynamic current is then proportional to the displacement of the WCs. Although the dynamic polarization is also gauge invariant only up to a “quantum of polarization,” the current J(t) is uniquely defined. In fact, it is the current J(t) which is the quantity that can be used for determining various properties of the system of interest as discussed in Sec. II E.

The theoretical framework laid out above is general and allows for arbitrarily strong and rapid variations of the homogeneous electric field. Consequently, one can perform RT-TDDFT calculations to study systems (both molecular and extended) under photoirradiation (i.e., photoexcitation), ion irradiation (i.e., electronic stopping), and static electric fields on an equal footing. As mentioned earlier, it should be noted that the underlying derivation of the dynamic polarization as a nonadiabatic Berry phase requires that the initial state be “Wannier-representable” (WR), as described in Ref. 65. Physically, Wannier-representability holds when the initial state of the system is insulating-like, not metallic, in the RT-TDDFT simulations. Additionally, in their work, Souza et al. proved that, in the absence of scattering, a WR state remains WR or “insulating-like” at all later times, even if the ground state of the Hamiltonian for a given ionic configuration becomes metallic.

With access to the dynamic current, one can obtain the frequency-dependent conductivity and also the dielectric function within linear response theory. For a system under a homogeneous perturbing field E(t) in the ν-direction, and for the time-dependent current in the μ-direction, the frequency dependent conductivity is obtained as

σ_{μ ν} (ω) = \frac{1}{\tilde{E_{ν}} (ω)} \int^{T} d t e^{i ω t} J_{μ} (t),

(18)

where $\tilde{E} (ω)$ is the Fourier transform of the applied electric field. For extended periodic systems, the frequency dependent dielectric function is related to the conductivity via

ε (ω) = 1 + \frac{4 π i}{3 ω} T r [σ_{μ ν} (ω)],

(19)

where Tr[σ(ω)] is the trace of the complex conductivity tensor. For extended systems, the imaginary part of this dielectric function is directly related to the optical absorption, whereas the real part is related to dispersion. For isolated systems, the macroscopic dielectric function is not well-defined, and instead the convention is to describe optical absorption in terms of the dipole strength function

S (ω) = \frac{4 π ω}{3 c} T r [I m σ_{μ ν} (ω)],

(20)

where σ_μν(ω) is generally referred to as the frequency-dependent polarizability in the case of isolated systems.⁹⁷ For computing absorption spectra using RT-TDDFT, one must choose an appropriate excitation procedure that simultaneously excites the system in a superposition of eigenstates. Any sudden perturbation at time t = 0 that is suddenly “switched off” at the next time step has the effect of inducing electronic oscillations in the system that includes all frequency components and thus results in a broadband electronic excitation of the system.⁹⁸ In RT-TDDFT simulations, there are two common choices for this sudden perturbation: first, there is the delta-function-like impulsive electric field. In this approach, a finite electric field is applied only for an infinitesimally small moment in time. In practice, however, RT-TDDFT calculations involve a finite time step Δt, which allows for an approximation of a true impulsive electric field. This has consequences for the electric field Fourier transform term $\tilde{E} (ω)$ ⁠, which for a true impulsive field with magnitude E₀ should yield $\tilde{E} (ω) = E_{0}$ ⁠. The equations defining the temporal profile and the Fourier transform of the electric field, with some magnitude E₀, are defined below for the impulse approach,

\begin{aligned} E_{i m p} (t) & = \{\begin{matrix} E_{0} & f o r t = 0 \\ 0 & f o r t > 0, \end{matrix}) \\ {\tilde{E}}_{i m p} (ω) & = E_{0} . \end{aligned}

(21)

In practice, however, due to the finite integration time step, the field profile is actually that of a “boxcar” function, and its Fourier transform results in a sinc function. An additional complication is that this approach can cause some difficulties in the numerical integration of the TDKS equations.⁹⁸ For these reasons, in this work, we primarily use a second, more numerically convenient “step function” electric field approach, proposed by Yabana et al.²³ In this alternative method, one performs a standard DFT calculation including a static uniform electric field to acquire a stationary state solution which is subsequently propagated in the TDKS equations for t > 0 without the field. This amounts to an adiabatic “switching on” of the field and a nonadiabatic “switching off” of the field at t = 0 of the RT-TDDFT simulation. The equations defining the temporal profile and the Fourier transform of the electric field, with some magnitude E₀, are defined below for the step function approach,

\begin{aligned} E_{s t e p} (t) & = \{\begin{matrix} E_{0} & f o r t < 0 \\ 0 & f o r t > 0, \end{matrix}) \\ {\tilde{E}}_{s t e p} (ω) & = \frac{E_{0}}{i ω} . \end{aligned}

(22)

In practice, of course, for a finite amount of time T. Consequently, the abrupt end of the oscillations in the polarization at the end of the numerical simulation leads to artificial “wiggles” in the calculated spectrum. The prominence of such numerical artifacts can be lessened to some degree by employing a damping function in the Fourier transform. A simple choice for the damping function, which we employ in this work, is the damping function f(t) = e^−γt, where γ is a damping constant. For RT-TDDFT simulations converged sufficiently with respect to total simulation time, the application of this damping function results in spectra with Lorentzian shaped peaks.

E. Implementation in QB@LL

In this work, we employ the plane-wave pseudopotential formalism in RT-TDDFT as discussed in Ref. 58. The simulations carried out in this work were performed using the highly parallelized plane-wave pseudopotential implementation of RT-TDDFT in the Qb@ll branch of the Qbox code,^58,59,99 in which we have implemented the real-time propagation of maximally localized Wannier functions. Details of the localization algorithm and its performance and accuracy are given in Subsections II E 1 and II E 2.

1. Maximal localization procedure

The method employed in this work for computing MLWFs and time-dependent MLWFs (TD-MLWFs) is an extension of the method developed by Gygi et al.¹⁰⁰ Gygi et al. developed a parallelized Cardoso-Souloumiac diagonalization algorithm for the efficient and accurate computation of MLWFs in molecular and extended systems. While their implementation was restricted to real wavefunctions (i.e., static ground state wavefunctions), the algorithm itself, as suggested by the authors, could be extended to complex wavefunctions (i.e., time-dependent wavefunctions), and in this work, we have carried out such an extension. What follows is a brief outline of the MLWF algorithm generalized for complex wavefunctions. Further details can be found in the studies by Gygi et al.¹⁰⁰

In the seminal work by Resta,⁸⁶ the quantum-mechanical position operator is defined for periodic systems and its respective spread functional. In a rectangular periodic cell of dimensions L_x, L_y, L_z, at the center of the Brillouin zone, it can be shown that the spread associated with the position operator is equal to using the spread functional associated with a set of six self-adjoint operators $\{{\hat{A}}^{(k)}\}, k = 1, \dots, 6$ ⁠, defined as follows. For an N-electron system, and the set of KS (or TDKS) orbitals $\{ϕ_{i}\}, i = 1, \dots, N,$

σ_{\{{\hat{A}}^{(k)}\}}^{2} (\{ϕ_{i}\}) \equiv \sum_{k = 1}^{6} \sum_{i} [⟨ϕ_{i} |{\hat{A}}^{2}| ϕ_{i}⟩ - {⟨ϕ_{i} |\hat{A}| ϕ_{i}⟩}^{2}],

\begin{aligned} {\hat{A}}^{(1)} & \equiv \cos \frac{2 π}{L_{x}} x \\ {\hat{A}}^{(2)} & \equiv \sin \frac{2 π}{L_{x}} x \\ {\hat{A}}^{(3)} & \equiv \cos \frac{2 π}{L_{y}} y \\ {\hat{A}}^{(4)} & \equiv \sin \frac{2 π}{L_{y}} y \\ {\hat{A}}^{(5)} & \equiv \cos \frac{2 π}{L_{z}} z \\ {\hat{A}}^{(6)} & \equiv \sin \frac{2 π}{L_{z}} z . \end{aligned}

(23)

Minimization of $σ_{\{{\hat{A}}^{(k)}\}}^{2} (\{ϕ_{i}\})$ is achieved by calculating a unitary transformation matrix U that simultaneously maximally diagonalizes the matrices ${\hat{A}}^{(k)}$ ⁠. This simultaneous diagonalization is achieved via the Cardoso-Souloumiac algorithm.¹⁰¹ The diagonalization is carried out iteratively until the spread functional $σ_{\{{\hat{A}}^{(k)}\}}^{2}$ converges to a minimum within a chosen tolerance. The resultant diagonal elements of the matrices ${\hat{A}}^{(k)}$ can be used to determine the positions of the MLWF centers (MLWFCs) and the spreads of the MLWFs. The iterative implementation of the Cardoso-Souloumiac algorithm for the computation of this transformation matrix U is further detailed in the work by Gygi et al.,¹⁰⁰ with the differences for the case of real matrices noted.

Because the calculation of this unitary MLWF transformation (i.e., unitary localization operator) is carried out at each RT-TDDFT step (sometimes for thousands of simulation steps), computational performance is of great importance. Like the Jacobi algorithm, the Cardoso-Souloumiac algorithm is inherently parallel. The current implementation makes use of message passing interface (MPI) and high-performance linear algebra libraries (ScaLAPACK, BLACS) in conjunction with a processor-data rotation scheme¹⁰² for the efficient parallelized computation of the MLWF transformation. Scaling and performance data for the TD-MLWF implementation compared with the standard TDKS implementation of RT-TDDFT are shown in Sec. II E 2.

2. Scaling and accuracy

The calculation and application of the TD-MLWF unitary transformation during the RT-TDDFT propagation adds additional computational cost in the simulations stemming from the calculation of the localization operator. Thus, it is important to examine the performance and accuracy of the RT-TDDFT/TD-MLWF implementation relative to the standard RT-TDDFT/TDKS approach. As test systems, we chose two representative cases: an isolated benzene molecule and crystalline silicon. All of these simulations used the LDA exchange-correlation functional and Hamann-Schluter-Chiang-Vanderbilt (HSCV) norm-conserving pseudopotentials.¹⁰³ The adiabatic approximation was used for the time dependence of the exchange correlation functional.^104,105 First, standard DFT calculations were performed to acquire the ground state wavefunctions. Then, in order to create a perturbed nonequilibrium initial condition for the RT-TDDFT simulations, the atoms were translated by +0.01 bohr in the x-direction at the start of the time-propagation. For each system, we performed two simulations: a control case in which the TDKS equations are propagated with TDKS Bloch orbitals, and a test case involving TD-MLWF propagation. For all TD-MLWF simulations, a convergence tolerance of 10⁻⁸ was used in the joint approximate diagonalization algorithm described in Sec. II E 1.

The scalability tests for the benzene case were performed using a 30 Ry plane-wave kinetic cutoff energy, a 50 × 50 × 50 bohr simulation cell, and a 0.1 a.u. time step for 100 simulation steps with the enforced time-reversal symmetry (ETRS) integrator.¹⁰⁶ The simulations were performed using MPI on nodes with 44 Intel Xeon processors each. Figure 2 shows the results of the performance test, and for this small case, the TD-MLWF propagation approach does not add any appreciable computational cost until we approach several hundred cores. On 352 cores, a single TD-MLWF simulation step is 1.41 times the cost of a single TDKS simulation step. However, we also observe that the performance of the standard TDKS simulation also decreases when we reach over 500 cores. For such a small system (30 electrons), this is not surprising that we reach this performance bottleneck at only a few hundred cores.

FIG. 2.

View large Download slide

Performance results, indicated in time per RT-TDDFT step, showing the scalability of the TD-MLWF propagation (blue) vs the standard TDKS propagation (red) in RT-TDDFT simulations for an isolated benzene molecule (left) and 512-atom crystalline silicon supercell (right).

For testing scalability with the silicon crystal, a 512-atom supercell was used, with a 60 Ry plane-wave kinetic cutoff energy, a 0.1 a.u. time step for 100 simulation steps with the ETRS integrator. The simulations were performed using MPI on an IBM BG/Q system with 16 MPI tasks per node. The results, shown in Fig. 3, show that in all cases studied the TD-MLWF propagation comes at a greater computational cost than TDKS propagation, as expected. However, scaling up to larger numbers of cores, the performance gap slightly decreases between the two methods: on 256 cores, each TD-MLWF step is 1.88 times the cost, and on 32 768 cores, each TD-MLWF step is 1.88 times the cost. The TD-MLWF propagation simulations in these test simulations used very strict convergence criteria for the calculation of the localization operator, and with relaxation of these criteria, the computational cost of the TD-MLWF propagation could be reduced even further while maintaining good accuracy.

FIG. 3.

View large Download slide

RT-TDDFT simulation accuracy, measured in terms of energy drift, for simulations with two different integration schemes (ETRS and FORK) and two different types of orbital propagations (TD-MLWF and TDKS) as a function of number of integrations per a.u. of time for an isolated benzene molecule (left) and 64-atom crystalline silicon supercell (right).

The scalability results in this work show that the calculation and propagation of TD-MLWFs incur an additional computational cost, but that the additional cost does not significantly impede practical simulations. That being said, we have not yet explored the possibilities for TD-MLWFs to be used to actually reduce the cost of RT-TDDFT simulations in certain contexts. To our knowledge, this work represents the first use of MLWF propagation in RT-TDDFT, but the spatially localized characteristics of MLWFs have been exploited in the context of ground-state DFT calculations in insulating systems to achieve order-N calculations of exact exchange¹⁰⁷ as well as order-N first-principles molecular dynamics simulations based on a “divide and conquer” scheme.¹⁰⁸ The efficient calculation of hybrid exchange could allow for RT-TDDFT simulations with hybrid exchange-correlation functionals. Additionally, the “divide and conquer” method described by Osei-Kuffour and Fattebert¹⁰⁸ could potentially be used in large-scale RT-TDDFT simulations to improve the performance by avoiding global communications and possibly by limiting the number of orbitals that need to be propagated at each time step. We plan to investigate these possibilities in the future.

In addition to testing the scalability of our TD-MLWF implementation, we performed simulations to test its accuracy. In principle, when one transforms Bloch orbitals into MLWFs, all physical observables should remain unchanged because the transformation is unitary, preserving the total electron density. However, this assumes that the transformation matrix is exactly unitary. Numerically, however, the joint diagonalization algorithm used to compute the matrix at each time step is approximate and requires us to choose a convergence tolerance for the spread functional (see Sec. II E 1). Thus, it is important to compare the TD-MLWF and standard TDKS (Bloch orbitals) RT-TDDFT simulations to ensure that we can achieve excellent agreement in physical observables such as total electronic energy. Again, we chose two representative systems as test cases: an isolated benzene molecule in a 50 × 50 × 50 bohr simulation cell as well as a 64-atom crystalline silicon supercell. In order to examine the accuracy of TD-MLWF propagation with different numerical integrators, we performed simulations with both the fourth-order Runge-Kutta (FORK) method^58,109 and the ETRS method.¹⁰⁶ As before, we initialize the RT-TDDFT simulations by shifting all atomic coordinates by +0.01 bohr in the x direction relative to the coordinates used in the ground state DFT calculation. From this nonequilibrium starting point, the TDKS equations are propagated for 100 a.u. of time in two different ways: a control case in which the TDKS equations are propagated via Bloch orbitals and a test case in which the TDKS equations are propagated in terms of the TD-MLWFs. With all atoms frozen, the total energy serves as a constant of motion, meaning that any drift in the energy can be attributed to numerical errors in the real-time propagation. For the TD-MLWF cases, a convergence tolerance of 10⁻⁸ was used for the joint approximate diagonalization algorithm.

The results (see Fig. 3) indicate that there is significant dependence on the numerical integrator used (FORK vs ETRS), with the ETRS integrator conserving the total energy in the system better than FORK for a given Δt, especially for the case of the benzene molecule. In addition, the ETRS propagator can use much larger time steps than the FORK (0.25 a.u. vs 0.15 a.u.) before becoming completely unstable. However, in all cases, there is no significant difference (<0.1%) between the RT-TDDFT simulations with the standard Bloch TDKS propagation and the TD-MLWF propagation. This shows that our TD-MLWF implementation is capable of maintaining high numerical accuracy in the computation and application of the unitary gauge transformation at every time step.

III. RESULTS AND DISCUSSION

A. Optical excitation

1. Benzene molecule

As an example application for the real-time propagation of MLWFs, we have performed RT-TDDFT simulations on a benzene molecule in a vacuum to calculate the optical absorption spectrum. The optical properties of gas-phase benzene have been well-studied both in experimental and theoretical studies. In particular, the optical absorption spectrum has been calculated via both RT-TDDFT^38,110 and the Liouville-Lanczos approach to linear response TDDFT (LR-TDDFT),⁵⁷ and in the case of gaseous benzene, the calculated spectra are in very good agreement with reported experimental results, even with the LDA exchange-correlation approximation.^111,112

In our RT-TDDFT approach, a static polarized state electronic structure was first obtained adiabatically via a DFT calculation with an electric field of magnitude 0.001 a.u. present using the MLWF/length-gauge formulation described in Sec. II. This electric field magnitude is sufficiently small to ensure a linear response and to ensure oscillations of the electronic dipole that are large enough compared to any numerical noise. The LDA exchange-correlation (XC) functional was used in all simulations. The basis was expanded in plane waves up to a kinetic cutoff energy of 30 Ry. A large cubic simulation cell of 100 × 100 × 100 bohr was used, with the benzene ring lying in the xy plane. Such a large simulation cell was used in order to avoid interactions between periodic images of the system which can lead to extra “ripples” in the absorption spectrum for high-energy excitations where ionization becomes important.²³ Another approach to mitigate such finite size effects is to employ a complex absorbing potential at the simulation cell boundaries,²³ which is a capability that we plan to implement in the future. The RT-TDDFT simulations were initialized from the static field-polarized MLWFs. The enforced time reversal symmetry (ETRS)¹⁰⁶ was used for the numerical propagation of the TDKS equations. A time step of 0.1 a.u. was used. At the initial RT-TDDFT time step, the electric field is “turned off,” and the system propagates for 500 a.u. for each simulation. This sudden “switching off” of the electric field introduces a phase to the wavefunctions, resulting in continuous oscillations of the electronic dipole. Throughout the RT-TDDFT simulations, nuclear motion can be neglected due to the relatively short time scales and the small perturbations; thus, the nuclear positions are held fixed in order to reduce computational expense.

As shown in Fig. 4, the RT-TDDFT calculated spectrum shows good agreement with the experimental data. Indeed, other studies have shown that the optical absorption spectrum of molecular benzene can be successfully determined with TDDFT calculations,³⁸ even with the simple LDA approximation to the exchange-correlation functional.

FIG. 4.

View large Download slide

Dipole strength function (i.e., electronic absorption spectrum) for gas phase benzene calculated via TD-MLWFs/RT-TDDFT in this work (purple) compared to experimental data (black).

In addition to calculating the photoabsorption spectrum of benzene, the TD-MLWF approach allows us to decompose the spectrum into contributions from different chemical moieties. For molecular systems, the maximal localization transformation results in Wannier functions that are localized on bonds, and/or lone pairs. Thus, in the case of benzene, the TD-MLWFs and their corresponding centers can easily be associated with either carbon-carbon (CC) or carbon-hydrogen (CH) bonds. By Fourier-transforming the dynamic polarization of just the CC bonds or just the CH bonds, we can acquire the chemically decomposed spectra shown in Fig. 5. Interestingly, there are portions of the C–H spectrum which are negative. While negative values for the total absorption would be unphysical, the decomposed spectrum allows for this possibility. The negative “absorption” indicates that, at ∼7 eV, the oscillations of the CH bonds are out of phase with those of the CC bonds. This destructive interference results in total dynamic polarization oscillations that are smaller in magnitude than the CC bond dipole oscillations.

FIG. 5.

View large Download slide

Dipole strength function for gas-phase benzene (purple) decomposed via TD-MLWF centers into carbon-carbon bonds (blue) and carbon-hydrogen bonds (red). The negative values observed in the CH bond spectrum indicate oscillations that are out of phase with the CC bond oscillations.

2. Water molecule

Due to the availability of experimental absorption spectrum data over a wide range (∼0 eV to ∼200 eV), gaseous water is another interesting test case. LR-TDDFT based on Casida’s equation⁸ is the most widely used approach for computing the absorption spectra of materials and molecules. However, Casida’s equation is solved iteratively in a basis of (occupied) × (virtual) dimensions,¹⁸ meaning that calculations involving more than a few low-lying excited states become computationally prohibitive. Thus, the calculation of the broadband spectrum of gaseous water is one problem that lends itself to RT-TDDFT. Additionally, RT-TDDFT methods in plane-wave bases have the capabilities to capture high-energy excitations. This is unlike atom-centered basis set RT-TDDFT methods, which are known to exhibit spurious high-energy artifacts unless an imaginary molecular orbital-based absorbing potential is used.¹¹³

The simulation details closely follow those described for the benzene molecule case. A static polarized state electronic structure was first determined adiabatically via a DFT calculation with an electric field of magnitude 0.001 a.u. present. The Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) exchange-correlation functional¹¹⁴ was used in all simulations. The basis was expanded in plane waves up to a kinetic cutoff energy of 40 Ry. A cubic simulation cell of 100 × 100 × 100 bohr was used. The enforced time reversal symmetry (ETRS)¹⁰⁶ was used for the numerical propagation of the TDKS equations. A time step of 0.05 a.u. was used, with a total propagation of 250 a.u. for each simulation. At the initial RT-TDDFT time step, the system is perturbed by suddenly switching off the static homogeneous electric field. This perturbation introduces a phase to the wavefunctions, resulting in continuous oscillations of the electronic dipole. This electric field magnitude is sufficiently small to ensure a linear response, to conserve energy throughout the numerical propagation, and to acquire oscillations of the electronic dipole that are large compared to any numerical noise. The nuclei are held in fixed positions throughout the RT-TDDFT simulations.

The RT-TDDFT calculated spectrum alongside the experimental spectrum¹¹⁵ is shown in Fig. 6 in two different energy ranges, with the absorption spectra in a lower energy range (0 eV–30 eV) displaying three well-defined peaks in addition to the rapid onset of broadband absorption. While the RT-TDDFT results show good qualitative agreement with the experimental spectrum, there is a red-shift of the peaks by ∼2 eV. It is possible that RT-TDDFT with hybrid XC functionals or meta-GGA functionals could yield spectra in better agreement with experiments, and we plan to explore this in a future study. We also see that the broadband spectrum (0 eV–150 eV) is in overall agreement with the experimental data. This showcases the ability of the abilities of the RT-TDDFT simulations to acquire the absorption spectrum over a wide energy range, including photoemissionlike excitations.

FIG. 6.

View large Download slide

Dipole strength function (i.e., electronic absorption spectrum) for gas phase water calculated via TD-MLWFs/RT-TDDFT in this work (blue) compared to the experimental spectrum (black). The inset more clearly shows the low-energy portions of the spectra.

Similar to the molecular benzene case, the TD-MLWFs allow us to decompose the spectrum for a more detailed chemical understanding of the photoabsorption trends. The spectrum, decomposed in terms of contributions from oxygen-hydrogen (OH) bond TD-MLWFs and lone pair (LP) TD-MLWFs, is shown in Fig. 7. For the OH bond spectrum, we again see a region of negative magnitude, implying destructive interference in the coupling between OH bonds and LPs. In addition, we observe that the LP spectrum does not contain sharp peaks as seen in the OH bond spectrum, which comprises the three well-defined peaks in the 0–15 eV range. However, with regard to the overall magnitude, the LP excitations dominate at low energies below ∼15 eV, whereas the OH bond absorption dominates in a higher energy range.

FIG. 7.

View large Download slide

Dipole strength function for gas-phase water (blue) decomposed via TD-MLWF centers into oxygen-hydrogen bonds (light blue) and lone pairs (red). The negative values observed in the OH bond spectrum indicate oscillations that are out of phase with the lone pair bond oscillations.

3. Benzene in liquid water

With the capability to calculate dynamic polarizability in both molecular and condensed matter systems, we can use the TD-MLWF approach to study complex systems such as molecules solvated in liquid water. Known as the solvatochromic effect, the optical absorption spectrum of solute molecules, such as organic dyes, varies depending on the solvent.¹¹⁶ Thus, it is often important to account for solvent effects when performing excited state calculations on molecules, materials, and nanostructures. The most common approach to represent the solvent environment in TDDFT calculations is through a polarizable continuum model (PCM),¹¹⁷ which is computationally inexpensive but may not accurately describe solvatochromic effects. In addition, mixed quantum mechanical-classical (QM/MM) methods with standard or polarizable force fields can be used, providing a balance between accuracy and computational cost.^118–120 However, in QM/MM methods, the partitioning QM region can be a somewhat arbitrary or ambiguous choice. In addition, by construction, in QM/MM approaches, the MM region cannot be photoexcited. An alternative, albeit computationally expensive, approach is to treat the entire system (solute and solvent) quantum mechanically, performing TDDFT calculations on the full electronic structure.¹²¹ With the TD-MLWF approach, we can propagate the full electronic system with the TDKS equations, with the additional benefit that the orbitals are localized, allowing for the calculation of dynamic polarization for any individual molecule, or for the liquid system as a whole.

In order to demonstrate this approach, we again chose benzene as a test case. However, in this system, we explicitly solvated the benzene molecule in a cubic simulation cell (25.04 a.u. × 25.04 a.u. × 25.04 a.u.) containing 73 water molecules (614 electrons in total, including benzene). In order to acquire the absorption spectrum for a solvated molecule, one should in principle perform TDDFT calculations on a random ensemble of solvent/solute structures.¹²² However, due to the considerable computational cost of performing a RT-TDDFT calculation on even one structure containing over 600 electrons, we restrict ourselves here to two representative configurations known to be important in influencing the electronic structure of benzene,¹²³ for the purposes of illustrating the TD-MLWF utility. A first-principles molecular dynamics simulation using the CP2K code⁴⁸ with the strongly constrained and appropriately normed (SCAN) XC functional¹²⁴ and the double zeta for valence plus polarization function (DZVP) basis was used to acquire the two MD snapshots to use for the RT-TDDFT calculations. In configuration 1, the hydrogen of a single water molecule is directed toward the center of the benzene ring (see Fig. 8 top left). In configuration 2, there is no such “special” benzene-water molecule interaction (see Fig. 8 top right). For both configurations, a standard DFT calculation was performed to acquire the electronic ground states. Then, the RT-TDDFT simulations with TD-MLWF propagation were performed. For all simulations, the PBE exchange-correlation functional was used. Additional calculation parameters were 0.05 a.u. time step, 250 a.u. total simulation time, 40 Ry plane-wave cutoff energy, and 0.001 a.u. electric field impulse in the x, y, and z directions (separately) at t = 0.

FIG. 8.

View large Download slide

Absorption spectra for benzene in the water system (black) decomposed via TD-MLWF centers into contributions from the benzene molecule (red) and the water molecules (blue). Two configurations from two MD snapshots are represented: one structure with a water molecule hydrogen pointed into the benzene molecule pi cloud (left) and another without any such coordination (right).

Figure 8 shows the calculated absorption spectra of the entire benzene and water system with the dynamic polarizability obtained from the RT-TDDFT simulations. Using the positions of the TD-MLWF centers to decompose the dynamic polarization response, we also calculated the absorption spectrum of both the benzene molecule and liquid water separately. The total absorption spectrum does not have particularly well-separated peaks, but with this spectrum decomposition, the contributions from each molecular species become clearer. For the benzene molecule, there is a peak at ∼7 eV, as was also seen in the vacuum environment. However, for configuration 2, there is an additional low-energy peak at ∼5 eV. These results illustrate the usefulness and flexibility of the TD-MLWF approach in uncovering molecular-level details of complex chemical phenomena such as the solvatochromic effect.

4. Crystalline silicon

As an example of the TD-MLWF approach for periodic systems, we have performed RT-TDDFT simulations on crystalline silicon due to the fact that this material has been studied by a range of RT-TDDFT codes, results which provide useful points of comparison. For all simulations, the LDA exchange correlation functional was used. A plane-wave kinetic energy cutoff of 40 Ry was used. Additional calculation parameters are as follows: 0.05 a.u. time step, 250 a.u. total simulation time, 40 Ry plane-wave cutoff energy, and 0.001 a.u. electric field impulse strength.

For periodic crystalline systems, one must take particular care to ensure convergence with respect to finite size effects. Our current TD-MLWF implementation does not yet include capabilities to incorporate multiple k-points; thus, in this work, we use only the Γ point and a silicon supercell in order to acquire accurate results for the bulk material. Previous RT-TDDFT calculations on the optical properties of silicon have shown that k-point grids of 16 × 16 × 16 can be required to ensure the full convergence of the absorption spectrum.²⁶ Because of this, converging the spectrum with respect to supercell size may seem like a computationally intractable problem, as it would require a simulation cell with tens of thousands of electrons. However, we can take advantage of the lattice symmetry to greatly reduce the cost of obtaining a converged absorption spectrum. Due to the symmetry of the diamond cubic crystalline silicon, the linear response of the system will be identical for impulses in the x, y, and z directions. Thus, the diagonal element polarizability tensor [Eq. (18)] will all be equal, implying that only one RT-TDDFT simulation is required to obtain the full linear response of the system. With regard to the dimensions of the silicon supercell, as shown in the work by Darrigan et al.,¹²⁵ one can achieve convergence with respect to cell dimensions simply by extending the supercell in the direction of the applied electric field. Thus, if we are perturbing the system with an impulsive electric field in the x direction, we can simply repeat the 8-atom cell in the x direction until we reach convergence. Figure S2 (supplementary material) shows the resultant spectra for a range of supercells, and we observe that a 256-atom supercell is required to converge the imaginary part of the dielectric function.

In Fig. 9, we compare our results to those acquired via a range of other TDDFT codes. The SALMON (scalable ab initio light-matter simulator for optics and nanoscience) software package⁴⁹ involves a RT-TDDFT implementation in which the time-dependent Kohn-Sham orbitals are discretized on a 3D real-space grid. In the SALMON code, the perturbing field is represented in the velocity gauge by a shift of the vector potential and the dielectric function is determined with respect to the induced electric current density.⁴⁹ The SIESTA code employs RT-TDDFT calculations within a linear combination of atomic orbital (LCAO) basis set framework, and the finite field calculations are also carried out in the velocity-gauge, with a vector potential representation of the perturbation.²⁶ Figure 9 includes results from linear-response TDDFT calculations in a plane-wave pseudopotential implementation.¹²⁶ Although linear-response TDDFT represents a different theoretical framework altogether, the results should agree with those acquired through RT-TDDFT simulations, and, as shown in Fig. 9, they do. Generally speaking, all of the spectra are in good agreement, with the slight differences likely being attributable to numerical differences in the various implementations and basis sets. Although TDDFT is known to yield poor silicon spectra for conventional XC approximations, this example helps to illustrate the flexibility that real-time propagation of TD-MLWFs provides, eliminating any need to change gauge representations between isolated and periodic systems.

FIG. 9.

View large Download slide

Imaginary part of the dielectric function for crystalline silicon calculated in this work via TD-MLWFs/RT-TDDFT using the QB@LL code (blue). This spectrum is compared to reported results from LR-TDDFT (red), and two velocity-gauge RT-TDDFT calculations using the atomic orbital-based SIESTA code (orange) and the real-space-grid SALMON code (pink).

B. Quantized charge transport

Due to their unique connection to the Berry phase, Wannier functions are increasingly being used in the context of topological materials.¹²⁷ As posited by Resta and Vanderbilt, considering the behavior of Wannier centers (WCs) under a cyclic adiabatic evolution of the Hamiltonian can lead to some useful insights.⁸⁹ At the end of such a cyclic evolution, the WCs must return to their initial positions plus an integer multiple of the lattice constant R because the initial and final Hamiltonians are the same. As illustrated in a 1D model schematic (adapted from Resta and Vanderbilt⁸⁹) shown in Fig. 10, one can consider two possible evolutions of the WCs that satisfy this condition. The route involving a shift of the WCs by a lattice vector R corresponds to the quantized topological pumping phenomenon first discussed by Thouless.¹²⁸ In recent years, different types of Thouless pumps have been demonstrated experimentally,^129,130 and Wannier functions have been used as a formal means to understand the mechanisms of these dynamic topological phenomena. Most theoretical studies and descriptions of topological pumping assume complete adiabaticity of the Hamiltonian evolution. It is important to note here that topological Thouless pumping is a phenomenon that manifests from the complex phase of the wavefunctions, which means that any real MLWF method based on time-independent DFT would fail to exhibit such a topological phenomenon. Our TD-MLWF/RT-TDDFT implementation, however, in principle allows us to study the quantized charge transport phenomenon, including nonadiabatic effects and driving forces, from first principles. We plan to carry out detailed studies in the future.

FIG. 10.

View large Download slide

Schematic (adapted from Resta and Vanderbilt⁸⁹) illustrating the possible trajectories for MLWF centers under an adiabatic cyclic evolution of the Hamiltonian. The MLWFCs can either (a) return to their initial sites or (b) “hop” to a new site one lattice vector away.

Here, we demonstrate the ability for TD-MLWFs, and in particular their centers, to describe quantized charge transport in the semiconducting polymer system of trans polyacetylene. We use a quasimonochromatic electric field pulse applied along the axis of the polymer chain as a cyclic nonadiabatic driving force for the electronic system. Applying such a pulse to the electronic ground state of the system induces either continuous oscillations of the TD-MLWFCs or tunneling in, depending on the magnitude of the applied electric field. In order to drive quantized charge transport in one direction, we prepare a field-polarized resonant state as the initial condition of the RT-TDDFT simulation by first adiabatically “switching on” a static electric field. This static electric field has the effect of “tilting” the periodic potential well, breaking the energy landscape symmetry in the x-direction. Then, the field-polarized stat is subjected to a quasimonochromatic pulse, which then drives the quantized charge transport. In this way, we can drive quantized charge transport for a certain resonant frequency.

The simulation details for the first-principles charge transport simulations are as follows: a 28-atom trans-polyacetylene chain was contained in a periodic simulation cell with dimensions of (34.28 × 20.00 × 20.00 bohr), with the polymer aligned along the x axis. All atoms were represented by HSCV pseudopotentials. A plane-wave cutoff energy of 20 Ry was used. The LDA XC functional was used for all calculations. For the RT-TDDFT simulations, a 0.1 a.u. time step was used with the ETRS integrator. Two initial electronic states for the RT-TDDFT simulations were calculated (1) the electronic ground state calculated via DFT and (2) the field-polarized state calculated adiabatically via DFT in the presence of a static electric field in the +x direction with a magnitude of 2.5 × 10⁻³ a.u. For the RT-TDDFT simulations, these systems are then subjected to a quasimonochromatic electric field pulse with a frequency of 2.8 eV (calculated resonant frequency of double-bond TD-MLWFs), a maximum field strength of 1.0 × 10⁻³ a.u., which is enveloped in a Gaussian with a full-width at half-maximum of 1.0 eV. In Fig. 11, the temporal profiles of the electric fields are shown, with the adiabatic (DFT) and nonadiabatic (RT-TDDFT) portions of the simulations delineated.

FIG. 11.

View large Download slide

Temporal profiles of the homogeneous electric field for the two simulation cases. The region to the left of the dashed line, $t \in (- \infty, 0)$ ⁠, represents the adiabatic switching-on of the electric field (blue) in the time-independent DFT calculation. The region to the right $t \in [0,500] a.u.$ of the dashed line represents the beginning of the RT-TDDFT simulations with the quasimonochromatic pulse with (blue) and without (green) the static electric field.

The spatial x-axis coordinates of the double-bond TD-MLWFCs with respect to time are shown in Fig. 12 for both simulation cases. Although the electrons, and consequently the TD-MLWFs, are all in principle indistinguishable, in Fig. 12, we have highlighted some of the plotted data as a visual guide to understand the observed phenomena. In the pulse-only case (Fig. 12 left), the dynamics are fairly simple, with the quasimonochromatic pulse exciting a resonance in the double-bond TD-MLWFs that causes repetitive Rabi oscillations, but no net charge transport in one direction. In the presence of the static electric field in addition to the quasimonochromatic pulse (Fig. 12 right), however, the dynamics are clearly different. At a point in time near the end of the pulse (∼475 a.u.), there is a sudden jump in the TD-MLWFC positions. This discontinuity, an average of 1.81 bohr displacement in the −x direction, corresponds to double-bond TD-MLWFs tunneling from one site on the polyacetylene chain to another. Figure 13 illustrates this pictorially through the perspective of the chemical structure and by showing snapshots of the polyacetylene TD-MLWFs from one RT-TDDFT time step to the next, only 0.1 a.u. later. For both simulations, all nuclei were held in fixed positions. Thus, even after the applied electric field is turned off at t = 500 a.u., the electronic system continues to propagate freely and cyclically as in a Rabi oscillation. Consequently, we observe that the quantized charge transport, initially driven by the pulse, continues repetitively, regularly, on average every 67.9 a.u. of time and with −1.81 bohr displacement of each TD-MLWFC. The discretized tunneling phenomenon is also reflected in the TD-MLWF spread. Figure 14 shows the total change in spread of the TD-MLWFs throughout the RT-TDDFT simulations. For the pulse-only case, the dynamics are as expected, with the pulse having the effect of cyclic oscillations in the spread (i.e., oscillatory localization and delocalization of the TD-MLWFs). However, in the pulse + static field case, we observe again jump discontinuities simultaneous with the TD-MLWFC position jumps. As the TD-MLWFs come close to tunneling, the spread increases rapidly, and then suddenly decrease after the tunneling occurred and the TD-MLWFs return to being localized between a different pair of carbon atoms.

FIG. 12.

View large Download slide

Positions of the double-bond TD-MLWFCs on the polyacetylene chain (commensurate with the x-axis) throughout the RT-TDDFT simulations for the “pulse only” (left) and the “pulse + static field” (right) cases. The positions of the TD-MLWFCs are plotted as a point at each step in the RT-TDDFT simulation to show the presence of the jump discontinuities (right).

FIG. 13.

View large Download slide

(Upper) Skeletal chemical structure for the trans-polyacetylene chain illustrating the double-bond TD-MLWFCs hopping occurring in the RT-TDDFT simulation. (Lower) Two consecutive RT-TDDFT simulation snapshots showing the carbon and hydrogen atoms (gray and white) and the TD-MLWFCs (yellow and dark blue) as spheres. Here, the TD-MLWFC spheres have been scaled and colored proportionally to the spreads of the associated TD-MLWFs.

FIG. 14.

View large Download slide

Total change in the TD-MLWF spread throughout the RT-TDDFT simulations for the “pulse only” case (green dots) and the “pulse + static field” case (blue dots). The inset shows the zoomed-in view of a small time region to illustrate the sudden decrease in spread at the time of TD-MLWF tunneling.

C. Electronic stopping excitations

One physical process that RT-TDDFT has been useful in simulating is electronic stopping. Electronic stopping is the phenomenon that occurs when fast-moving charged particles (> ∼10 keV) transfer energy to the surrounding electronic system through excitation and ionization. Electronic stopping power is defined, as the rate of energy transfer from projectile ions to electrons is a key quantity in describing this process. In recent years, RT-TDDFT simulations have been quite successful in predicting the electronic stopping power of a range of matter and for a variety of projectile ions.^{16,30,32,35,36} However, there are more scientific insights to be gained from understanding ion irradiation excitation dynamics on an atomistic level. Induced density analysis¹⁶ and single-particle projection methods²⁷ have been previously used as tools to understand the molecular level details of the stopping process. However, both of these approaches have limitations. Analysis in terms of the time-dependent induced density can be difficult to interpret, and quantitative comparisons usually require the use of a charge-partitioning scheme, which relies on a reasonable, but still rather arbitrary, choice. Projection methods, on the other hand, also have their limitations. In particular, with projections onto the equilibrium (i.e., ground state) eigenstates, it is difficult to include enough virtual/unoccupied eigenstates in the projection calculation to account for excited electrons in ionization, which is rather common in the electronic stopping process. Thus, in practice, projection methods are only most useful in analyzing excited hole distributions. With these challenges in mind, TD-MLWF propagation provides a new, useful tool. By construction, TD-MLWFs can give molecular level details of the excited electrons in the electronic stopping process without the need to resort to the charge density partitioning scheme. Additionally, a key challenge in studies of electronic stopping is to examine differences in the electron excitation induced by photon and ion radiation. The TD-MLWF approach, with its flexibility in simulations of both types of radiations, can aid in this challenge.

In order to demonstrate some of the capabilities of TD-MLWFs for descriptions of electronic stopping processes, we have performed several RT-TDDFT simulations of proton irradiation of an isolated benzene molecule. The simulations were performed with a 0.1 a.u. time step in the ETRS propagator, a 30 Ry plane-wave kinetic energy cutoff, the LDA XC functional, and HSCV pseudopotentials. The setup was as follows: a proton with constant velocity in the +x direction was moved through a 30 × 30 × 50 bohr simulation cell containing a benzene molecule oriented in the xy plane. Four proton velocities were simulated (v = 0.625, 1.25, 5.00, 10.0 a.u.) at two different impact parameters relative to the plane of the benzene molecule. Once the proton has traversed the length of the simulation cell (50 bohr), it is removed, at which point the RT-TDDFT simulation and TD-MLWF propagation are continued for 500 a.u. The electronic dipole moment of the excited benzene molecule, calculated via the TD-MLWFCs, then continues to oscillate. Taking the absolute value of the Fourier transform of the electronic dipole moment presents a spectrum of the excitations caused in the proton electronic stopping process. Figure 15 shows the resultant “electronic stopping spectra,” alongside the calculated optical spectrum. For the low-velocity cases (0.625 and 1.25 a.u.), the spectra primarily comprise a single well-defined peak at ∼6.9 eV, which aligns with the lowest energy peak in the optical spectrum, indicating that at lower velocities, the proton primarily excites the mode that corresponds to the first optical excitation. Approaching higher proton velocities (5.00 and 10.0 a.u.), the peak at ∼6.9 eV diminishes, and we noted the emerging peaks in higher energy ranges (15 eV–40 eV). This result is consistent with previous RT-TDDFT studies of electronic stopping which have shown that more ionization occurs as the ion velocity increases. We also observed impact parameter dependence in the electronic stopping spectra (see Fig. 15, left panel compared to right panel). Not surprisingly, there is an overall trend of increasing spectrum intensity with decreasing impact parameter for all velocities. However, the low-velocity spectra (0.625 a.u. and 1.25 a.u.) show large increases in intensity (by a factor or ∼2), whereas the high-velocity spectra show much smaller increases.

FIG. 15.

View large Download slide

“Electronic stopping spectra” for a proton impinging on a gas phase benzene molecule with an impact parameter of 10 a.u. (left) and an impact parameter of 8 a.u. (right). Four velocities (0.625, 1.25, 5.00, and 10.0 a.u.) were simulated, and four spectra were calculated for each case. The optical spectrum (dashed line) is shown for comparison.

IV. CONCLUSIONS

Real-time, time-dependent density functional theory (RT-TDDFT) has gained popularity as a first-principles approach to study a variety of excited-state phenomena such as optical excitations and electronic stopping. Within RT-TDDFT simulations, the gauge freedom of the time-dependent electronic orbitals can be exploited for numerical and scientific convenience while the unitary transformation does not alter physical properties calculated from the quantum dynamics of electrons. Exploiting this gauge freedom, we demonstrated the propagation of maximally localized Wannier functions (MLWFs) in RT-TDDFT within our plane-wave pseudopotential RT-TDDFT branch of the QB@LL code.

In recent years, the gauge freedom in RT-TDDFT has been exploited in various ways including using the so-called parallel transport gauge for numerical efficiency.⁶² The time-dependent MLWF (TD-MLWF) formalism detailed in this work represents another such gauge representation that is practical for a wide range of RT-TDDFT simulations. Due to the MLWF connection to the Berry phase, the dynamic polarization of both isolated molecular systems and condensed matter systems in periodic boundary conditions is well-defined, allowing for a scalar potential representation of homogeneous electric fields. This length gauge formulation circumvents the commonly used velocity gauge²⁶ formulation (within the electrodynamics) in which the electric field needs to be represented as a magnetic flux when periodic boundary conditions are used. As a demonstration, we have performed RT-TDDFT simulations of the linear response of both molecular and crystalline systems in order to calculate the resultant absorption spectra. For the case of benzene, we observed good agreement between the RT-TDDFT calculated spectrum and experimental absorption spectrum. As another molecular case, we used the TD-MLWF method to calculate the absorption spectrum of an isolated (i.e., gas phase) water molecule. In this case, we observed excellent qualitative agreement with the experimental optical absorption spectrum covering a wide range of excitation energies (0–150 eV). Notably, the simulations do not give rise to any spurious peaks in the high energy region, a boon that is likely a result of the plane-wave pseudopotential basis upon which the RT-TDDFT code is implemented. Additionally, we have demonstrated the utility of TD-MLWFs for analyzing intermolecular and intramolecular excitation dynamics. In the case of the isolated molecules (benzene and water), we demonstrated the decomposition of the optical absorption spectra in terms of different chemical moieties, such as lone-pair oxygen-hydrogen bonds and carbon-carbon bonds. We showed how this TD-MLWF decomposition scheme can also be used to examine the absorption spectrum of a benzene molecule solvated in liquid water, showing the possibilities for investigating solvatochromic effects. As a capability demonstration of the length-gauge formulation for representing the homogeneous electric field in periodic boundary conditions, we calculated the dielectric function of crystalline silicon, a system which has been well-studied with other TDDFT codes. Our spectrum shows a good agreement with both LR-TDDFT results¹²⁶ and RT-TDDFT results from a variety of codes in atomic orbital²⁶ and real-space grid⁴⁹ implementations.

One application of RT-TDDFT that was not discussed in this work is the computation of core-level optical spectra, such as X-ray absorption spectra (XAS). While we have recently studied the core electron excitation dynamics in electronic stopping successfully,¹³¹ there remain significant challenges in the computation of XAS with RT-TDDFT for a number of reasons. In the context of atomic-orbital based RT-TDDFT simulations, such calculations have been carried out in recent years.^24,25,41,132 However, with plane-wave pseudopotential (PW-PP) approaches, the accurate representation of core electron wavefunctions requires an extremely high plane-wave cutoff (i.e., PW basis set expansion) which makes the calculation of XAS particularly computationally difficult in practice. In general, for any underlying basis set of choice, there remains shortcoming with regard to the exchange-correlation approximation, with popular hybrid functionals even failing to accurately predict XAS excitation energies.¹³² Although the TD-MLWF approach can also be used for modeling of XAS using RT-TDDFT, it is a topic that requires a thorough investigation in a future work.

We presented a few examples of the utility of the TD-MLWF gauge representation beyond the linear response regime (e.g., optical absorption spectrum). Using the ability to simulate both static and time-dependent electric fields as scalar potentials in the length gauge, we demonstrated the simulation of quantized charge transport in polyacetylene. Although it is a relative simple test case, similar simulations could be used in the future to carry first-principles investigations into dynamic topological transport phenomena, beyond the adiabatic dynamic description of an ideal topological Thouless pump.^128,133 Finally, we presented TD-MLWF RT-TDDFT simulations of the electronic stopping of a proton impinging on a benzene molecule. Although it is often difficult to make direct comparisons between the excitation behavior of systems under photon and ion irradiation, we have proposed a type of analysis through the calculation of the “electronic stopping spectrum” based on the dynamic polarization; in the same way, the optical absorption spectrum is calculated from the dynamics polarization. Using the electronic stopping spectrum, one can identify what parts of the optical absorption spectrum are excited in an electronic stopping process as a function of the projectile proton velocity.

There are many additional applications of RT-TDDFT where TD-MLWFs could provide advantages over conventional TDKS propagation. In the context of ground-state DFT, the spatial locality of MLWFs has been used to greatly accelerate calculations with hybrid exchange-correlation (XC) functionals¹⁰⁷ and first principles molecular dynamics simulations.¹⁰⁸ We have demonstrated that, for large systems, the TD-MLWF implementation scales well over tens of thousands of cores without incurring prohibitive computational costs (i.e., only two times more expensive). Thus, the TD-MLWF approach makes it possible to also accelerate hybrid XC calculations in the context of RT-TDDFT for studying large complex condensed matter systems.

SUPPLEMENTARY MATERIAL

See supplementary material for information on vacuum regions and absorbing potentials for molecular optical spectra, silicon supercell convergence, and quantized charge transport pulse frequency dependence.

ACKNOWLEDGMENTS

The authors would like to thank Erik W. Draeger for helpful discussions. This work was supported by the National Science Foundation under Grants Nos. CHE-1565714, DGE-1144081, and OAC-17402204. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357.

REFERENCES

1.

P.

Hohenberg

and

W.

Kohn

,

Phys. Rev.

136

,

B864

(

1964

).

https://doi.org/10.1103/physrev.136.b864

Google Scholar

Crossref

2.

M.

Grätzel

,

J. Photochem. Photobiol., C

4

,

145

(

2003

).

https://doi.org/10.1016/s1389-5567(03)00026-1

Google Scholar

Crossref

3.

A. J.

Nozik

,

Annu. Rev. Phys. Chem.

52

,

193

(

2001

).

https://doi.org/10.1146/annurev.physchem.52.1.193

Google Scholar

Crossref

PubMed

4.

O. I.

Obolensky

,

Nucl. Instrum. Methods Phys. Res. B

266

,

1623

(

2008

).

https://doi.org/10.1016/j.nimb.2007.11.054

Google Scholar

Crossref

5.

E.

Runge

and

E. K.

Gross

,

Phys. Rev. Lett.

52

,

997

(

1984

).

https://doi.org/10.1103/physrevlett.52.997

Google Scholar

Crossref

6.

M. A.

Marques

and

E. K.

Gross

,

Annu. Rev. Phys. Chem.

55

,

427

(

2004

).

https://doi.org/10.1146/annurev.physchem.55.091602.094449

Google Scholar

Crossref

PubMed

7.

C. A.

Ullrich

,

Time-Dependent Density-Functional Theory: Concepts and Applications

(

Oxford University Press

,

Oxford

,

2011

).

Google Scholar

Crossref

8.

M. E.

Casida

,

Recent Advances in Density Functional Methods (Part I)

(

World Scientific

,

1995

), pp.

155

–

192

.

Google Scholar

Crossref

9.

R. L.

Martin

,

J. Chem. Phys.

118

,

4775

(

2003

).

https://doi.org/10.1063/1.1558471

Google Scholar

Crossref

10.

M.

Petersilka

,

U. J.

Gossmann

, and

E. K. U.

Gross

,

Phys. Rev. Lett.

76

,

1212

(

1996

).

https://doi.org/10.1103/physrevlett.76.1212

Google Scholar

Crossref

PubMed

11.

K.

Burke

,

J.

Werschnik

, and

E. K. U.

Gross

,

J. Chem. Phys.

123

,

062206

(

2005

).

https://doi.org/10.1063/1.1904586

Google Scholar

Crossref

12.

C. A.

Ullrich

and

A. D.

Bandrauk

, in

Time-Dependent Density Functional Theory

, Lecture Notes in Physics, edited by

M. A. L.

Marques

,

N. T.

Maitra

,

F. M. S.

Nogueira

,

E. K. U.

Gross

, and

A.

Rubio

(

Springer-Verlag

,

Berlin

,

2012

), Vol. 837, pp.

351

–

371

.

Google Scholar

Crossref

13.

V. A.

Goncharov

,

J. Chem. Phys.

139

,

084104

(

2013

).

https://doi.org/10.1063/1.4818807

Google Scholar

Crossref

PubMed

14.

X.

Andrade

,

S.

Botti

,

M. A.

Marques

, and

A.

Rubio

,

J. Chem. Phys.

126

,

184106

(

2007

).

https://doi.org/10.1063/1.2733666

Google Scholar

Crossref

PubMed

15.

F.

Ding

,

B. E. V.

Kuiken

,

B. E.

Eichinger

, and

X.

Li

,

J. Chem. Phys.

138

,

064104

(

2013

).

https://doi.org/10.1063/1.4790583

Google Scholar

Crossref

PubMed

16.

D. C.

Yost

and

Y.

Kanai

,

Phys. Rev. B

94

,

115107

(

2016

).

https://doi.org/10.1103/physrevb.94.115107

Google Scholar

Crossref

17.

J.

Liu

and

J. M.

Herbert

,

J. Chem. Phys.

143

,

034106

(

2015

).

https://doi.org/10.1063/1.4926837

Google Scholar

Crossref

PubMed

18.

R. E.

Stratmann

,

G. E.

Scuseria

, and

M. J.

Frisch

,

J. Chem. Phys.

109

,

8218

(

1998

).

https://doi.org/10.1063/1.477483

Google Scholar

Crossref

19.

S.

Baroni

and

R.

Gebauer

, in

Fundamentals of Time-Dependent Density Functional Theory

, Lecture Notes in Physics, edited by

M.

Marques

,

N.

Maitra

,

F.

Nogueira

,

E.

Gross

, and

A.

Rubio

(

Springer

,

Berlin, Heidelberg

,

2012

), Vol. 837, pp.

375

–

390

.

Google Scholar

Crossref

20.

W.

Liang

,

S. A.

Fischer

,

M. J.

Frisch

, and

X.

Li

,

J. Chem. Theory Comput.

7

,

3540

(

2011

).

https://doi.org/10.1021/ct200485x

Google Scholar

Crossref

PubMed

21.

J. J.

Goings

,

P. J.

Lestrange

, and

X.

Li

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

8

,

e1341

(

2018

).

https://doi.org/10.1002/wcms.1341

Google Scholar

Crossref

22.

K.

Yabana

and

G. F.

Bertsch

,

Phys. Rev. B

54

,

4484

(

1996

).

https://doi.org/10.1103/physrevb.54.4484

Google Scholar

Crossref

23.

K.

Yabana

,

T.

Nakatsukasa

,

J. I.

Iwata

, and

G. F.

Bertsch

,

Phys. Status Solidi B

243

,

1121

(

2006

).

https://doi.org/10.1002/pssb.200642005

Google Scholar

Crossref

24.

K.

Lopata

,

B. E. V.

Kuiken

,

M.

Khalil

, and

N.

Govind

,

J. Chem. Theory Comput.

8

,

3284

(

2012

).

https://doi.org/10.1021/ct3005613

Google Scholar

Crossref

PubMed

25.

R. G.

Fernando

,

M. C.

Balhoff

, and

K.

Lopata

,

J. Chem. Theory Comput.

11

,

646

(

2015

).

https://doi.org/10.1021/ct500943m

Google Scholar

Crossref

PubMed

26.

C. D.

Pemmaraju

,

F. D.

Vila

,

J. J.

Kas

,

S. A.

Sato

,

J. J.

Rehr

,

K.

Yabana

, and

D.

Prendergast

,

Comput. Phys. Commun.

226

,

30

(

2018

).

https://doi.org/10.1016/j.cpc.2018.01.013

Google Scholar

Crossref

27.

K. G.

Reeves

and

Y.

Kanai

,

Sci. Rep.

7

,

40379

(

2017

).

https://doi.org/10.1038/srep40379

Google Scholar

Crossref

PubMed

28.

Y.

Miyamoto

and

A.

Rubio

,

J. Phys. Soc. Jpn.

87

,

041016

(

2018

).

https://doi.org/10.7566/jpsj.87.041016

Google Scholar

Crossref

29.

S.

Meng

and

E.

Kaxiras

,

Nano Lett.

10

,

1238

(

2010

).

https://doi.org/10.1021/nl100442e

Google Scholar

Crossref

PubMed

30.

A.

Schleife

,

Y.

Kanai

, and

A. A.

Correa

,

Phys. Rev. B

91

,

014306

(

2015

).

https://doi.org/10.1103/physrevb.91.014306

Google Scholar

Crossref

31.

A.

Ojanperä

,

A. V.

Krasheninnikov

, and

M.

Puska

,

Phys. Rev. B

89

,

035120

(

2014

).

https://doi.org/10.1103/physrevb.89.035120

Google Scholar

Crossref

32.

A. A.

Correa

,

Comput. Mater. Sci.

150

,

291

(

2018

).

https://doi.org/10.1016/j.commatsci.2018.03.064

Google Scholar

Crossref

33.

J. M.

Pruneda

,

D.

Sánchez-Portal

,

A.

Arnau

,

J. I.

Juaristi

, and

E.

Artacho

,

Phys. Rev. Lett.

99

,

235501

(

2007

).

https://doi.org/10.1103/physrevlett.99.235501

Google Scholar

Crossref

PubMed

34.

I.

Maliyov

,

J.-P.

Crocombette

, and

F.

Bruneval

,

Eur. Phys. J. B

91

,

172

(

2018

).

https://doi.org/10.1140/epjb/e2018-90289-y

Google Scholar

Crossref

35.

D. C.

Yost

,

Y.

Yao

, and

Y.

Kanai

,

Phys. Rev. B

96

,

115134

(

2017

).

https://doi.org/10.1103/physrevb.96.115134

Google Scholar

Crossref

36.

K. G.

Reeves

,

Y.

Yao

, and

Y.

Kanai

,

Phys. Rev. B

94

,

041108

(

2016

).

https://doi.org/10.1103/physrevb.94.041108

Google Scholar

Crossref

37.

K.

Lopata

and

N.

Govind

,

J. Chem. Theory Comput.

7

,

1344

(

2011

).

https://doi.org/10.1021/ct200137z

Google Scholar

Crossref

PubMed

38.

K.

Yabana

and

G. F.

Bertsch

,

Int. J. Quantum Chem.

75

,

55

(

1999

).

https://doi.org/10.1002/(sici)1097-461x(1999)75:1<55::aid-qua6>3.0.co;2-k

Google Scholar

Crossref

39.

M. A.

Marques

,

X.

López

,

D.

Varsano

,

A.

Castro

, and

A.

Rubio

,

Phys. Rev. Lett.

90

,

258101

(

2003

).

https://doi.org/10.1103/physrevlett.90.258101

Google Scholar

Crossref

PubMed

40.

A. R.

Attar

,

A.

Bhattacherjee

,

C. D.

Pemmaraju

,

K.

Schnorr

,

K. D.

Closser

,

D.

Prendergast

, and

S. R.

Leone

,

Science

356

,

54

(

2017

).

https://doi.org/10.1126/science.aaj2198

Google Scholar

Crossref

PubMed

41.

C. D.

Pemmaraju

,

Comput. Condens. Matter

18

,

e00348

(

2018

).

https://doi.org/10.1016/j.cocom.2018.e00348

Google Scholar

Crossref

42.

J. J.

Goings

and

X.

Li

,

J. Chem. Phys.

144

,

234102

(

2016

).

https://doi.org/10.1063/1.4953668

Google Scholar

Crossref

PubMed

43.

B.

Peng

,

D. B.

Lingerfelt

,

F.

Ding

,

C. M.

Aikens

, and

X.

Li

,

J. Phys. Chem. C

119

,

6421

(

2015

).

https://doi.org/10.1021/acs.jpcc.5b00263

Google Scholar

Crossref

44.

C. M.

Isborn

,

X.

Li

, and

J. C.

Tully

,

J. Chem. Phys.

126

,

134307

(

2007

).

https://doi.org/10.1063/1.2713391

Google Scholar

Crossref

PubMed

45.

C. L.

Moss

,

C. M.

Isborn

, and

X.

Li

,

Phys. Rev. A

80

,

024503

(

2009

).

https://doi.org/10.1103/physreva.80.024503

Google Scholar

Crossref

46.

Y.

Miyamoto

,

H.

Zhang

,

X.

Cheng

, and

A.

Rubio

,

Phys. Rev. B

96

,

115451

(

2017

).

https://doi.org/10.1103/physrevb.96.115451

Google Scholar

Crossref

47.

Y.

Takimoto

,

F. D.

Vila

, and

J. J.

Rehr

,

J. Chem. Phys.

127

,

154114

(

2007

).

https://doi.org/10.1063/1.2790014

Google Scholar

Crossref

PubMed

48.

J.

Hutter

,

M.

Iannuzzi

,

F.

Schiffmann

, and

J.

VandeVondele

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

4

,

15

(

2014

).

https://doi.org/10.1002/wcms.1159

Google Scholar

Crossref

49.

M.

Noda

,

S. A.

Sato

,

Y.

Hirokawa

,

M.

Uemoto

,

T.

Takeuchi

,

S.

Yamada

,

A.

Yamada

,

Y.

Shinohara

,

M.

Yamaguchi

,

K.

Iida

et al.

Comput. Phys. Commun.

235

,

356

(

2019

).

https://doi.org/10.1016/j.cpc.2018.09.018

Google Scholar

Crossref

50.

X.

Andrade

,

J.

Alberdi-Rodriguez

,

D. A.

Strubbe

,

M. J.

Oliveira

,

F.

Nogueira

,

A.

Castro

,

J.

Muguerza

,

A.

Arruabarrena

,

S. G.

Louie

,

A.

Aspuru-Guzik

et al.

J. Phys.: Condens. Matter

24

,

233202

(

2012

).

https://doi.org/10.1088/0953-8984/24/23/233202

Google Scholar

Crossref

PubMed

51.

A.

Castro

,

H.

Appel

,

M.

Oliveira

,

C. A.

Rozzi

,

X.

Andrade

,

F.

Lorenzen

,

M. A.

Marques

,

E. K. U.

Gross

, and

A.

Rubio

,

Phys. Status Solidi B

243

,

2465

(

2006

).

https://doi.org/10.1002/pssb.200642067

Google Scholar

Crossref

52.

T. S.

Nguyen

and

J.

Parkhill

,

J. Chem. Theory Comput.

11

,

2918

(

2015

).

https://doi.org/10.1021/acs.jctc.5b00262

Google Scholar

Crossref

PubMed

53.

W.

Liang

,

C. T.

Chapman

, and

X.

Li

,

J. Chem. Phys.

134

,

184102

(

2011

).

https://doi.org/10.1063/1.3589144

Google Scholar

Crossref

PubMed

54.

X.

Li

,

S. M.

Smith

,

A. N.

Markevitch

,

D. A.

Romanov

,

R. J.

Levis

, and

H. B.

Schlegel

,

Phys. Chem. Chem. Phys.

7

,

233

(

2005

).

https://doi.org/10.1039/b415849k

Google Scholar

Crossref

PubMed

55.

M. J.

Frisch

,

G. W.

Trucks

,

H. B.

Schlegel

,

G. E.

Scuseria

,

M. A.

Robb

,

J. R.

Cheeseman

,

G.

Scalmani

,

V.

Barone

,

G. A.

Petersson

,

H.

Nakatsuji

et al. gaussian 16, Revision B.01,

Gaussian, Inc.

,

Wallingford, CT

,

2016

.

56.

F.

Bruneval

,

T.

Rangel

,

S. M.

Hamed

,

M.

Shao

,

C.

Yang

, and

J. B.

Neaton

,

Comput. Phys. Commun.

208

,

149

(

2016

).

https://doi.org/10.1016/j.cpc.2016.06.019

Google Scholar

Crossref

57.

B.

Walker

,

A. M.

Saitta

,

R.

Gebauer

, and

S.

Baroni

,

Phys. Rev. Lett.

96

,

113001

(

2006

).

https://doi.org/10.1103/physrevlett.96.113001

Google Scholar

Crossref

PubMed

58.

A.

Schleife

,

E. W.

Draeger

,

Y.

Kanai

, and

A. A.

Correa

,

J. Chem. Phys.

137

,

22A546

(

2012

).

https://doi.org/10.1063/1.4758792

Google Scholar

Crossref

PubMed

59.

F.

Gygi

,

IBM J. Res. Dev.

52

,

137

(

2008

).

https://doi.org/10.1147/rd.521.0137

Google Scholar

Crossref

60.

A.

Schleife

,

E. W.

Draeger

,

V. M.

Anisimov

,

A. A.

Correa

, and

Y.

Kanai

,

Comput. Sci. Eng.

16

,

54

(

2014

).

https://doi.org/10.1109/mcse.2014.55

Google Scholar

Crossref

61.

G. F.

Bertsch

,

J. I.

Iwata

,

A.

Rubio

, and

K.

Yabana

,

Phys. Rev. B

62

,

7998

(

2000

).

https://doi.org/10.1103/physrevb.62.7998

Google Scholar

Crossref

62.

W.

Jia

,

D.

An

,

L.-W.

Wang

, and

L.

Lin

,

J. Chem. Theory Comput.

14

,

5645

(

2018

).

https://doi.org/10.1021/acs.jctc.8b00580

Google Scholar

Crossref

PubMed

63.

D.

An

and

L.

Lin

, preprint arXiv:1804.02095 (

2018

).

64.

W.

Jia

and

L.

Lin

,

Comput. Phys.

(

2019

).

https://doi.org/10.1016/j.cpc.2019.02.009

Google Scholar

65.

I.

Souza

,

J.

Íniguez

, and

D.

Vanderbilt

,

Phys. Rev. B

69

,

085106

(

2004

).

https://doi.org/10.1103/physrevb.69.085106

Google Scholar

Crossref

66.

E. K.

Gross

,

J. F.

Dobson

, and

M.

Petersilka

, in

Density Functional Theory II: Topics in Current Chemistry

, edited by

R. F.

Nalewajski

(

Springer

,

Berlin, Heidelberg

,

1996

), Vol. 181, pp.

81

–

172

.

Google Scholar

Crossref

67.

J. I.

Fuks

,

L.

Lacombe

,

S. E.

Nielsen

, and

N. T.

Maitra

,

Phys. Chem. Chem. Phys.

20

,

26145

(

2018

).

https://doi.org/10.1039/c8cp03957g

Google Scholar

Crossref

PubMed

68.

G.

Vignale

,

C. A.

Ullrich

, and

S.

Conti

,

Phys. Rev. Lett.

79

,

4878

(

1997

).

https://doi.org/10.1103/physrevlett.79.4878

Google Scholar

Crossref

69.

C. A.

Ullrich

and

K.

Burke

,

J. Chem. Phys.

121

,

28

(

2004

).

https://doi.org/10.1063/1.1756865

Google Scholar

Crossref

PubMed

70.

C. A.

Ullrich

,

J. Chem. Phys.

125

,

234108

(

2006

).

https://doi.org/10.1063/1.2406069

Google Scholar

Crossref

PubMed

71.

H. O.

Wijewardane

and

C. A.

Ullrich

,

Phys. Rev. Lett.

100

,

056404

(

2008

).

https://doi.org/10.1103/physrevlett.100.056404

Google Scholar

Crossref

PubMed

72.

C.-W.

Lee

and

A.

Schleife

,

Eur. Phys. J. B

91

,

222

(

2018

).

https://doi.org/10.1140/epjb/e2018-90204-8

Google Scholar

Crossref

73.

R.

Ullah

,

E.

Artacho

, and

A. A.

Correa

,

Phys. Rev. Lett.

121

,

116401

(

2018

).

https://doi.org/10.1103/physrevlett.121.116401

Google Scholar

Crossref

PubMed

74.

F.

Bloch

,

Z. Phys.

52

,

555

(

1929

).

https://doi.org/10.1007/bf01339455

Google Scholar

Crossref

75.

G. H.

Wannier

,

Phys. Rev.

52

,

191

(

1937

).

https://doi.org/10.1103/physrev.52.191

Google Scholar

Crossref

76.

W.

Kohn

,

Phys. Rev.

115

,

809

(

1959

).

https://doi.org/10.1103/physrev.115.809

Google Scholar

Crossref

77.

J. D.

Cloizeaux

,

Phys. Rev.

129

,

554

(

1963

).

https://doi.org/10.1103/physrev.129.554

Google Scholar

Crossref

78.

N.

Marzari

,

A. A.

Mostofi

,

J. R.

Yates

,

I.

Souza

, and

D.

Vanderbilt

,

Rev. Mod. Phys.

84

,

1419

(

2012

).

https://doi.org/10.1103/revmodphys.84.1419

Google Scholar

Crossref

79.

U.

Stephan

,

R. M.

Martin

, and

D. A.

Drabold

,

Phys. Rev. B

62

,

6885

(

2000

).

https://doi.org/10.1103/physrevb.62.6885

Google Scholar

Crossref

80.

W.

Ku

,

H.

Rosner

,

W. E.

Pickett

, and

R. T.

Scalettar

,

Phys. Rev. Lett.

89

,

167204

(

2002

).

https://doi.org/10.1103/physrevlett.89.167204

Google Scholar

Crossref

PubMed

81.

W.-C.

Lu

,

C.-Z.

Wang

,

T.-L.

Chan

,

K.

Ruedenberg

, and

K.-M.

Ho

,

Phys. Rev. B

70

,

041101

(

2004

).

https://doi.org/10.1103/physrevb.70.041101

Google Scholar

Crossref

82.

V. P.

Smirnov

and

D. E.

Usvyat

,

Phys. Rev. B

64

,

245108

(

2001

).

https://doi.org/10.1103/physrevb.64.245108

Google Scholar

Crossref

83.

B.

Sporkmann

and

H.

Bross

,

Phys. Rev. B

49

,

10869

(

1994

).

https://doi.org/10.1103/physrevb.49.10869

Google Scholar

Crossref

84.

B.

Sporkmann

and

H.

Bross

,

J. Phys.: Condens. Matter

9

,

5593

(

1997

).

https://doi.org/10.1088/0953-8984/9/26/008

Google Scholar

Crossref

85.

N.

Marzari

and

D.

Vanderbilt

,

Phys. Rev. B

56

,

12847

(

1997

).

https://doi.org/10.1103/physrevb.56.12847

Google Scholar

Crossref

86.

R.

Resta

,

Phys. Rev. Lett.

80

,

1800

(

1998

).

https://doi.org/10.1103/physrevlett.80.1800

Google Scholar

Crossref

87.

R.

Resta

,

Rev. Mod. Phys.

66

,

899

(

1994

).

https://doi.org/10.1103/revmodphys.66.899

Google Scholar

Crossref

88.

R. D.

King-Smith

and

D.

Vanderbilt

,

Phys. Rev. B

47

,

1651

(

1993

).

https://doi.org/10.1103/physrevb.47.1651

Google Scholar

Crossref

89.

R.

Resta

and

D.

Vanderbilt

,

Physics of Ferroelectrics

, Topics in Applied Physics (

Springer

,

Berlin, Heidelberg

,

2007

), Vol. 105, pp.

31

–

68

.

Google Scholar

Crossref

90.

I.

Souza

,

J.

Íniguez

, and

D.

Vanderbilt

,

Phys. Rev. Lett.

89

,

117602

(

2002

).

https://doi.org/10.1103/physrevlett.89.117602

Google Scholar

Crossref

PubMed

91.

R. W.

Nunes

and

X.

Gonze

,

Phys. Rev. B

63

,

155107

(

2001

).

https://doi.org/10.1103/physrevb.63.155107

Google Scholar

Crossref

92.

N. T.

Maitra

,

I.

Souza

, and

K.

Burke

,

Phys. Rev. B

68

,

045109

(

2003

).

https://doi.org/10.1103/physrevb.68.045109

Google Scholar

Crossref

93.

K.

Yabana

,

T.

Sugiyama

,

Y.

Shinohara

,

T.

Otobe

, and

G. F.

Bertsch

,

Phys. Rev. B

85

,

045134

(

2012

).

https://doi.org/10.1103/physrevb.85.045134

Google Scholar

Crossref

94.

M.

Stengel

and

N. A.

Spaldin

,

Phys. Rev. B

75

,

205121

(

2007

).

https://doi.org/10.1103/physrevb.75.205121

Google Scholar

Crossref

95.

M.

Sharma

,

Y.

Wu

, and

R.

Car

,

Int. J. Quantum Chem.

95

,

821

(

2003

).

https://doi.org/10.1002/qua.10633

Google Scholar

Crossref

96.

D.

Vanderbilt

and

R. D.

King-Smith

,

Phys. Rev. B

48

,

4442

(

1993

).

https://doi.org/10.1103/physrevb.48.4442

Google Scholar

Crossref

97.

A.

Zangwill

and

P.

Soven

,

Phys. Rev. A

21

,

1561

(

1980

).

https://doi.org/10.1103/physreva.21.1561

Google Scholar

Crossref

98.

C. A.

Ullrich

,

Time-Dependent Density-Functional Theory: Concepts and Applications

(

Oxford University Press

,

Great Clarendon Street, Oxford

,

2012

), pp.

204

–

210

.

Google Scholar

99.

E. W.

Draeger

,

X.

Andrade

,

J. A.

Gunnels

,

A.

Bhatele

,

A.

Schleife

, and

A. A.

Correa

,

J. Parallel Distrib. Comput.

106

,

205

(

2017

).

https://doi.org/10.1016/j.jpdc.2017.02.005

Google Scholar

Crossref

100.

F.

Gygi

,

J.-L.

Fattebert

, and

E.

Schwegler

,

Comput. Phys. Commun.

155

,

1

(

2003

).

https://doi.org/10.1016/s0010-4655(03)00315-1

Google Scholar

Crossref

101.

J.-F.

Cardoso

and

A.

Souloumiac

,

SIAM J. Matrix Anal. Appl.

17

,

161

(

1996

).

https://doi.org/10.1137/s0895479893259546

Google Scholar

Crossref

102.

G. H.

Golub

and

C. F. V.

Loan

, in

Matrix Computations

(

JHU Press

,

1996

).

Google Scholar

103.

D.

Vanderbilt

,

Phys. Rev. B

32

,

8412

(

1985

).

https://doi.org/10.1103/physrevb.32.8412

Google Scholar

Crossref

104.

A.

Zangwill

and

P.

Soven

,

Phys. Rev. Lett.

45

,

204

(

1980

).

https://doi.org/10.1103/physrevlett.45.204

Google Scholar

Crossref

105.

A.

Zangwill

and

P.

Soven

,

Phys. Rev. B

24

,

4121

(

1981

).

https://doi.org/10.1103/physrevb.24.4121

Google Scholar

Crossref

106.

A.

Castro

,

M. A. L.

Marques

, and

A.

Rubio

,

J. Chem. Phys.

121

,

3425

(

2004

).

https://doi.org/10.1063/1.1774980

Google Scholar

Crossref

PubMed

107.

X.

Wu

,

A.

Selloni

, and

R.

Car

,

Phys. Rev. B

79

,

085102

(

2009

).

https://doi.org/10.1103/physrevb.79.085102

Google Scholar

Crossref

108.

D.

Osei-Kuffuor

and

J.-L.

Fattebert

,

Phys. Rev. Lett.

112

,

046401

(

2014

).

https://doi.org/10.1103/physrevlett.112.046401

Google Scholar

Crossref

PubMed

109.

W. H.

Press

,

S. A.

Teukolsky

,

W. T.

Vetterling

, and

B. P.

Flannery

,

Numerical Recipes: The Art of Scientific Computing

, 3rd ed. (

Cambridge University Press

,

Cambridge

,

2007

).

Google Scholar

110.

X.

Qian

,

J.

Li

,

X.

Lin

, and

S.

Yip

,

Phys. Rev. B

73

,

035408

(

2006

).

https://doi.org/10.1103/physrevb.73.035408

Google Scholar

Crossref

111.

E.-E.

Koch

and

A.

Otto

,

Chem. Phys. Lett.

12

,

476

(

1972

).

https://doi.org/10.1016/0009-2614(72)90011-5

Google Scholar

Crossref

112.

A.

Hiraya

and

K.

Shobatake

,

J. Chem. Phys.

94

,

7700

(

1991

).

https://doi.org/10.1063/1.460155

Google Scholar

Crossref

113.

K.

Lopata

and

N.

Govind

,

J. Chem. Theory Comput.

9

,

4939

(

2013

).

https://doi.org/10.1021/ct400569s

Google Scholar

Crossref

PubMed

114.

J. P.

Perdew

,

K.

Burke

, and

M.

Ernzerhof

,

Phys. Rev. Lett.

77

,

3865

(

1996

).

https://doi.org/10.1103/physrevlett.77.3865

Google Scholar

Crossref

PubMed

115.

H.

Hayashi

,

N.

Watanabe

,

Y.

Udagawa

, and

C. C.

Kao

,

J. Chem. Phys.

108

,

823

(

1998

).

https://doi.org/10.1063/1.475446

Google Scholar

Crossref

116.

A.

Marini

,

A.

Munoz-Losa

,

A.

Biancardi

, and

B.

Mennucci

,

J. Phys. Chem. B

114

,

17128

(

2010

).

https://doi.org/10.1021/jp1097487

Google Scholar

Crossref

PubMed

117.

M.

Cossi

,

V.

Barone

,

R.

Cammi

, and

J.

Tomasi

,

Chem. Phys. Lett.

255

,

327

(

1996

).

https://doi.org/10.1016/0009-2614(96)00349-1

Google Scholar

Crossref

118.

J. M.

Olsen

,

K.

Aidas

,

K. V.

Mikkelsen

, and

J.

Kongsted

,

J. Chem. Theory Comput.

6

,

249

(

2009

).

https://doi.org/10.1021/ct900502s

Google Scholar

Crossref

119.

K.

Sneskov

,

T.

Schwabe

,

O.

Christiansen

, and

J.

Kongsted

,

Phys. Chem. Chem. Phys.

13

,

18551

(

2011

).

https://doi.org/10.1039/c1cp22067e

Google Scholar

Crossref

PubMed

120.

C. M.

Isborn

,

A. W.

Gotz

,

M. A.

Clark

,

R. C.

Walker

, and

T. J.

Martínez

,

J. Chem. Theory Comput.

8

,

5092

(

2012

).

https://doi.org/10.1021/ct3006826

Google Scholar

Crossref

PubMed

121.

J. M.

Milanese

,

M. R.

Provorse

,

E.

Alameda

, Jr., and

C. M.

Isborn

,

J. Chem. Theory Comput.

13

,

2159

(

2017

).

https://doi.org/10.1021/acs.jctc.7b00159

Google Scholar

Crossref

PubMed

122.

T. J.

Zuehlsdorff

and

C. M.

Isborn

,

Int. J. Quantum Chem.

119

,

e25719

(

2019

).

https://doi.org/10.1002/qua.25719

Google Scholar

Crossref

123.

M.

Allesch

,

F. C.

Lightstone

,

E.

Schwegler

, and

G.

Galli

,

J. Chem. Phys.

128

,

014501

(

2008

).

https://doi.org/10.1063/1.2806288

Google Scholar

Crossref

PubMed

124.

J.

Sun

,

A.

Ruzsinszky

, and

J. P.

Perdew

,

Phys. Rev. Lett.

115

,

036402

(

2015

).

https://doi.org/10.1103/physrevlett.115.036402

Google Scholar

Crossref

PubMed

125.

C.

Darrigan

,

M.

Rérat

,

G.

Mallia

, and

R.

Dovesi

,

J. Comput. Chem.

24

,

1305

(

2003

).

https://doi.org/10.1002/jcc.10274

Google Scholar

Crossref

PubMed

126.

F.

Sottile

,

F.

Bruneval

,

A. G.

Marinopoulos

,

L. K.

Dash

,

S.

Botti

,

V.

Olevano

,

N.

Vast

,

A.

Rubio

, and

L.

Reining

,

Int. J. Quantum Chem.

102

,

684

(

2005

).

https://doi.org/10.1002/qua.20486

Google Scholar

Crossref

127.

D.

Gresch

,

G.

Autes

,

O. V.

Yazyev

,

M.

Troyer

,

D.

Vanderbilt

,

B. A.

Bernevig

, and

A. A.

Soluyanov

,

Phys. Rev. B

95

,

075146

(

2017

).

https://doi.org/10.1103/physrevb.95.075146

Google Scholar

Crossref

128.

D. J.

Thouless

,

Phys. Rev. B

27

,

6083

(

1983

).

https://doi.org/10.1103/physrevb.27.6083

Google Scholar

Crossref

129.

S.

Nakajima

,

T.

Tomita

,

S.

Taie

,

T.

Ichinose

,

H.

Ozawa

,

L.

Wang

,

M.

Troyer

, and

Y.

Takahashi

,

Nat. Phys.

12

,

296

(

2016

).

https://doi.org/10.1038/nphys3622

Google Scholar

Crossref

130.

M.

Lohse

,

C.

Schweizer

,

O.

Zilberberg

,

M.

Aidelsburger

, and

I.

Bloch

,

Nat. Phys.

12

,

350

(

2016

).

https://doi.org/10.1038/nphys3584

Google Scholar

Crossref

131.

Y.

Yao

,

D. C.

Yost

, and

Y.

Kanai

, preprint arXiv:1903.03194 (

2019

).

132.

P. J.

Lestrange

,

P. D.

Nguyen

, and

X.

Li

,

J. Chem. Theory Comput.

11

,

2994

(

2015

).

https://doi.org/10.1021/acs.jctc.5b00169

Google Scholar

Crossref

PubMed

133.

L.

Privitera

,

A.

Russomanno

,

R.

Citro

, and

G. E.

Santoro

,

Phys. Rev. Lett.

120

,

106601

(

2018

).

https://doi.org/10.1103/physrevlett.120.106601

Google Scholar

Crossref

PubMed

2019

Author(s)

Propagation of maximally localized Wannier functions in real-time TDDFT

I. INTRODUCTION

II. THEORETICAL AND COMPUTATIONAL DETAILS

A. Time-dependent Kohn Sham equations

B. Bloch and Wannier orbital representations

C. Finite electric fields and periodic boundary conditions

D. Dynamic polarization via Wannier functions

E. Implementation in QB@LL

1. Maximal localization procedure

2. Scaling and accuracy

III. RESULTS AND DISCUSSION

A. Optical excitation

1. Benzene molecule

2. Water molecule

3. Benzene in liquid water

4. Crystalline silicon

B. Quantized charge transport

C. Electronic stopping excitations

IV. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Propagation of maximally localized Wannier functions in real-time TDDFT

I. INTRODUCTION

II. THEORETICAL AND COMPUTATIONAL DETAILS

A. Time-dependent Kohn Sham equations

B. Bloch and Wannier orbital representations

C. Finite electric fields and periodic boundary conditions

D. Dynamic polarization via Wannier functions

E. Implementation in QB@LL

1. Maximal localization procedure

2. Scaling and accuracy

III. RESULTS AND DISCUSSION

A. Optical excitation

1. Benzene molecule

2. Water molecule

3. Benzene in liquid water

4. Crystalline silicon

B. Quantized charge transport

C. Electronic stopping excitations

IV. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only