Imaginary-time time-dependent density functional theory and its application for robust convergence of electronic states

TL;DR

This paper introduces an imaginary-time evolution method within density functional theory to achieve more reliable and robust convergence to electronic ground states, especially in challenging systems where traditional SCF methods struggle.

Contribution

The paper presents a novel imaginary-time TDDFT approach that guarantees convergence to the ground state regardless of the exchange-correlation functional used, improving robustness over conventional methods.

Findings

Successfully applied to systems with convergence issues

Maintains self-consistency throughout the process

Theoretically justified via the van Leeuwen theorem

Abstract

Reliable and robust convergence to the electronic ground state within density functional theory (DFT) Kohn-Sham (KS) calculations remains a thorny issue in many systems of interest. In such cases, charge sloshing can delay or completely hinder the convergence. Here, we use an approach based on transforming the time-dependent DFT equations to imaginary time, followed by imaginary-time evolution, as a reliable alternative to the self-consistent field (SCF) procedure for determining the KS ground state. We discuss the theoretical and technical aspects of this approach and show that the KS ground state should be expected to be the long-imaginary-time output of the evolution, independent of the exchange-correlation functional or the level of theory used to simulate the system. By maintaining self-consistency between the single-particle wavefunctions and the electronic density throughout the determination of the stationary state, our method avoids the typical difficulties encountered in SCF. To demonstrate dependability of our approach, we apply it to selected systems which struggle to converge with SCF schemes. In addition, through the van Leeuwen theorem, we affirm the physical meaningfulness of imaginary time TDDFT, justifying its use in certain topics of statistical mechanics such as in computing imaginary time path integrals.

Equations

$$|{\Phi }_i(t)⟩=e^{-\text{i}t\hat{H}}|{\Phi }_i⟩=e^{-\text{i}t{E}_i}|{\Phi }_i⟩,$$

$$|{\Phi }_i(\tau )⟩=e^{-\tau \hat{H}}|{\Phi }_i⟩=e^{-\tau E_i}|{\Phi }_i⟩.$$

$$|\ast ⟩\tilde{\Psi }(\tau )=\sum _{i=0}^{\infty }{A}_i(0)e^{-\tau E_i}|{\Phi }_i⟩,$$

$$|\Psi (\tau )⟩=\sum _{i=0}^{\infty }\frac{A_i(0)e^{-\tau E_i}}{\Omega (\tau )}|{\Phi }_i⟩,$$

$$\text{i}\text{partialderivative}t{\varphi }_j(\mathbf{r},t)$$

$${\hat{H}}_{\text{KS}}[n(\mathbf{r},t)]$$

$$v_s[n(\mathbf{r},t)]=v(\mathbf{r})+v_{\text{H}}\text{quantity}[n(\mathbf{r},t)]+v_{\text{xc}}\text{quantity}[n\text{quantity}(\mathbf{r},t)].$$

$$v_{\text{H}}\text{quantity}[n(\mathbf{r},t)]$$

$$v_{\text{xc}}\text{quantity}[n\text{quantity}(\mathbf{r},t)]$$

$$n(\mathbf{r},t)$$

$$|{\varphi }_j(t)⟩$$

$$\hat{U}\text{quantity}(t,t_0)$$

$$|{\varphi }_j(\tau )⟩$$

$$\hat{\mathcal{U}}\text{quantity}(\tau ,{\tau }_0)$$

$$\hat{\mathcal{U}}\text{quantity}(\tau +\Delta \tau ,\tau )$$

$${\hat{H}}_{\text{KS}}(\tau )$$

$$\hat{H}(t)=\hat{T}+\hat{V}(t)+\hat{W}+\lambda (t),$$

$$\hat{T}$$

$$\hat{V}(t)$$

$$\hat{W}$$

$$|\Psi (t)⟩$$

$$\hat{U}\text{quantity}(t,t_0)$$

$$={\hat{U}}_{\lambda }(t,t_0)\hat{\mathcal{T}}{\hat{U}}_{\hat{H}}(t,t_0),$$

$${\hat{U}}_{\lambda }(t,t_0)$$

$${\hat{U}}_{\hat{H}}(t,t_0)$$

$$\hat{\mathcal{U}}\text{quantity}(\tau ,{\tau }_0)$$

$${\hat{\mathcal{U}}}_{\lambda }(\tau ,{\tau }_0)$$

$${\hat{\mathcal{U}}}_{\hat{H}}(\tau ,{\tau }_0)$$

$${\hat{\mathcal{U}}}_{\lambda }(\tau ,{\tau }_0)=\text{quantity}[\sum _{j=0}^{\infty }\text{absolutevalue}{A_j(0)}^2{e}^{-2\tau E_j}{]}^{-1/2}$$

$$\lambda (\tau )=\frac{1}{2}\text{derivative}\tau \ln \text{quantity}[\sum _{j=0}^{\infty }\text{absolutevalue}{A_j(0)}^2{e}^{-2\tau E_j}]=-\text{expectationvalue}E(\tau ),$$

$$\hat{\mathcal{U}}\text{quantity}(\tau ,{\tau }_0)$$

$${\hat{\mathcal{U}}}_{{\lambda }_{\text{KS}}}(\tau ,{\tau }_0)$$

$$D[n,n_0]\equiv \frac{1}{2}{d}_1(n,n_0)=\frac{1}{2}\int \text{absolutevalue}n(\mathbf{r})-n_0(\mathbf{r})\text{differential}[3]\mathbf{r},$$

Full text

I Introduction

Density functional theory (DFT) is a widely used approach enabling ab initio calculations of electronic and material properties. Unlike direct approaches to studying quantum systems through the Schrödinger equation where the wavefunction is the central object, DFT uses the electron density $n(\mathbf{r})$ as the fundamental physical quantity. In principle, through the Hohenberg-Kohn theoremHohenberg and Kohn (1964), the ground state $n(\mathbf{r})$ of a system uniquely determines all of its observables. It is standard practice to use the Kohn-Sham (KS) systemKohn and Sham (1965) of non-interacting fermions as a shortcut to obtaining the ground state density, employing specially formulated potentialsCeperley and Alder (1980); Perdew et al. (1996) that are functionals of $n(\mathbf{r})$ to approximate the electron-electron interactions.

There are many techniques to find the ground state of the KS equations, including methods which: (a) aim at direct determination of the minimum of the KS total energy functionalŠtich et al. (1989); Voorhis and Head-Gordon (2002); Weber et al. (2008); and (b) use iterative methods based on diagonalization of the KS Hamiltonian in conjunction with iterative improvements of the ground state charge density through mixing. The present work is distinct from both these approaches. We focus on comparing our method to type (b) approaches.

For a system with $N$ electrons, the lowest $N$ eigenstates to the KS equations determine $n(\mathbf{r})$, which itself appears in the KS equations through an effective single-particle potential. In general, finding the set of $N$ eigenstates that satisfy the KS equations involves an iterative process known as self-consistent field (SCF) iterations that produce successively better approximations to the solution. In its simplest conceptualization the iterative approach involves solving the eigenvalue problem for an initial density distribution, then using the resulting eigenstates to produce the next approximation to the density. When this approach is iterated, except for the simplest systems, it rarely converges to a self-consistent solution. In order to stabilize the SCF loops and improve the convergence rate, various mixing schemes are typically employed. These schemes take advantage of the information contained in multiple previous trial densities to select the next one. A popular mixing scheme is direct inversion of the iterative subspace (DIIS), also known as Pulay mixingPulay (1980, 1982).

When SCF schemes require many iterations to reach an acceptable solution, or fail to converge, the choices are to change the mixing scheme or its parameters, start with a different density, or fractionally occupy statesMichelini et al. (1998); Rabuck and Scuseria (1999) which some methods implement by introducing a fictitious electronic temperature (Fermi smearingMermin (1965); Verstraete and Gonze (2001)). If these fail, one can resort to computationally-intensive direct minimization methodsŠtich et al. (1989); Voorhis and Head-Gordon (2002); Weber et al. (2008) to find a solution. The convergence difficulties for SCF usually arise in systems with large unit cells and in metallic systemsAnglade and Gonze (2008), or when an excited state is desired. The small differences in eigenenergies of the single-particle states, as well as the presence of many states near the Fermi level, can cause very different eigenstates to be occupied from step to step. This can lead to large variations in the density, causing the phenomenon known as charge sloshingKresse and Furthmüller (1996) where a fluctuating charge density from step to step is observed with insufficient attenuation to reach convergence.

In the present paper we transform the time-dependent KS (TDKS) equations of time-dependent density functional theory (TDDFT)Peuckert (1978); Runge and Gross (1984); van Leeuwen (2001) to imaginary timeKolesov et al. (2018). We use these equations to propagate an initial state to very long imaginary time, refining it down to the KS state corresponding to its lowest energy component. The idea of using imaginary-time propagation (ITP) to find eigenstates is well-known, and it is frequently used to find ground state solutions to the Schrödinger equation describing single-particle systems with a fixed potentialBader et al. (2013); Lehtovaara et al. (2007). Imaginary time steps have also been used to find self-consistent solutions to the Hartree-Fock equationsDavies et al. (1980) and for nuclear energy density functional calculationsRyssens et al. (2015). It has also been employed in a DFT context as an alternative to the diagonalization step to find the single-particle KS eigenstates for a fixed electronic densityHernández et al. (2007); Aichinger and Krotscheck (2005). However, imaginary-time evolution has yet to be examined as a stand-alone substitute to iterative density updating in solving the KS equations. In the present method both the density and wavefunction evolve together towards the ground state according to the imaginary time TDKS equations, remaining consistent with each other throughout the calculation. We discuss the theoretical foundation of the imaginary-time evolution of the KS system, a procedure which is non-unitary, requiring re-orthonormalization of the states at each imaginary time-step. We show that the proof provided by van Leeuwenvan Leeuwen (2001) for TDDFT can be extended to imaginary-time TDDFT (it-TDDFT), affirming in principle that the density of a KS system will evolve in imaginary time in the same manner as the true many-body interacting system. The imaginary-time propagation method in DFT has attractive theoretical and practical benefits when applied to systems that are challenging to study using standard methods of solving the KS equations, as we demonstrate on model systems.

We benchmark our approach by applying it to the benzene molecule and show that it converges to the same ground state energy as other SCF-based methods. Next, we apply our method to systems with known difficulties in achieving convergence. We chose to examine a copper nanocluster Cu13 with fixed magnetization and a spin-unpolarized Ru55 nanocluster. We show that self-consistent solutions are hard to realize in both systems using the most popular standard approach, SCF with Pulay mixing. In general, we find that while requiring more computation, our method is more dependable and more autonomous compared to SCF. It provides a good alternative to existing methodologies when the latter fail to converge in challenging systems, or if a user wishes to find an unfamiliar system’s ground state with minimal intervention; this can be particularly useful when computations are carried out in an automated fashion on large clusters of processors.

The paper is organized as follows: Section II presents the method, Section III extends the van Leeuwen theorem to imaginary time and discusses certain theoretical considerations that establish the method’s robustness, Section IV gives some example calculations, and Section V contains our conclusions.

II Methodology

II.1 Imaginary-Time Propagation

First, let us take the Hamiltonian $\hat{H}$ to be time-independent. Under the substitution $t\to-\textrm{i}\tau$, where $\tau$ is real, the time evolution operator transforms from $e^{-\textrm{i}t\hat{H}}$ to $e^{-\tau\hat{H}}$. When $\ket{\Phi_{i}}$ is an eigenstate of $\hat{H}$, the time evolution of the eigenstate switches from rotating its phase proportionally to its energy:

[TABLE]

to shrinking its amplitude by an exponential factor with a rate proportional to the energy:

[TABLE]

For the case of a time-dependent Hamiltonian the previous equations still hold for infinitesimal amounts of time $\Delta t$ or $\Delta\tau$. For an arbitrary initial wavefunction $\ket{\Psi(0)}$, imaginary-time propagation amounts to

[TABLE]

where $A_{i}(0)$ is the amplitude of the eigenstate component initially present. As imaginary time goes to infinity, $\tau\to\infty$, the eigenstate $\ket{\Phi_{j}}$ corresponding to the lowest energy eigenvalue with $A_{j}(0)\neq 0$ will dominate. We can choose to keep the state $\ket{\Psi(\tau)}$ normalized by dividing by the norm $\Omega(\tau)\equiv\sqrt{\innerproduct*{\tilde{\Psi}(\tau)}{\tilde{\Psi}(\tau)}}=\sqrt{\sum_{i=0}^{\infty}\absolutevalue{A_{i}(0)}^{2}e^{-2\tau E_{i}}}$,

[TABLE]

which then yields $\lim_{\tau\to\infty}\ket{\Psi(\tau)}=\ket{\Phi_{j}}$. Note that the state $\ket{\Psi(\tau)}$ could refer to a single-particle wavefunction $\Psi(\mathbf{r})$ or a many-body wavefunction $\Psi(\mathbf{r}_{1},\ldots,{\mathbf{r}}_{N})$; the above discussion is applicable to either case. Since an arbitrary initial state generated by randomizing the coefficients in some basis is likely to have a nonzero ground state component, ITP is often used to find ground state wavefunctions and energies.

II.2 Implementation within the Kohn-Sham Formalism

In TDDFT, starting from an initial state, the KS system obeys the equations of motion (in atomic units):

[TABLE]

with time-dependent effective potential

[TABLE]

In these expressions, $v({\mathbf{r}})$ is the external potential and

[TABLE]

The Kohn-Sham time-evolution can be reformulated in terms of a time-propagator which acts on single-particle states and is given by

[TABLE]

where $\hat{\mathcal{T}}$ is the time-ordering operator. In imaginary time, applying the substitution $t\to-\textrm{i}\tau$ results in

[TABLE]

where $\hat{\mathcal{T}}_{\tau}$ now time-orders in imaginary time. Note that the imaginary-time propagator is not unitary.

Employing the same numerical scheme used for real time propagation of KS states on an atomic basisKolesov et al. (2016), we evolve in imaginary-time the single-particle states using finite time steps $\Delta\tau$ and we approximate the instantaneous imaginary-time propagator with the second-order Magnus expansion:

[TABLE]

The Hamiltonian at the midpoint is approximated as the average of the Hamiltonians at $\tau_{i}$ and $\tau_{i+1}$, $\hat{H}_{\textrm{KS}}\quantity(\tau_{i}+\frac{\Delta\tau}{2})\approx\frac{1}{2}\quantity[\hat{H}_{\textrm{KS}}(\tau_{i})+\hat{H}_{\textrm{KS}}(\tau_{i+1})]$. Each step is iterated to self-consistency in order to make use of the Hamiltonian at $\tau_{i+1}$. We use the Padé rational polynomial approximation of arbitrary degree to obtain the general matrix exponential. Further details of the numerical propagation can be found in our earlier workKolesov et al. (2016), which describes TDAP-2.0, a TDDFT code we used, built on top of SIESTASoler et al. (2002), a DFT package which uses strictly localized basis sets. While the midpoint Hamiltonian greatly aids stability and energy conservation in real time propagation, in practice we have found that for imaginary-time propagation we can just use the first step in the iterative procedure, which simply applies the approximation $\hat{H}_{\textrm{KS}}\quantity(\tau_{i}+\frac{\Delta\tau}{2})\approx\hat{H}_{\textrm{KS}}\quantity(\tau_{i})$. This explicit propagation is faster since the Hamiltonian only needs to be evaluated once per propagation step, and the effect on the size of the maximum stable time-step appears negligible compared to the implicit method using the midpoint Hamiltonian. This is expected since imaginary-time propagation is inherently more stable than the real-time propagation the TDAP-2.0 code was originally designed to solve.

Because the imaginary-time propagator is not unitary, the single-particle states lose their normalization and generally cease to be orthogonal. The simple expression for density in Eq. (9) becomes more complicated if the single-particle states $\phi_{j}$ are non-orthonormal. It is convenient to reorthonormalize the single-particle states at each time step. The details of how the orthogonalization is achieved do not affect the physics, as we show in Section III.2. We use the modified Gram-Schmidt algorithm to orthonormalize the states.

While we employ a localized atomic basis for our calculations, the method we propose is independent of the basis used to represent the Kohn-Sham orbitals, and can easily be implemented in other popular bases, like plane waves or gaussians.

III Theoretical Considerations

III.1 Van Leeuwen Theorem in Imaginary Time

The van Leeuwen theorem states that a time-dependent particle density $n({\mathbf{r}},t)$ belonging to a many-particle system with two-particle interaction $\hat{W}$ can always be reproduced by a unique (up to an additive purely time-dependent constant) external potential $v^{\prime}({\mathbf{r}},t)$ in another many-particle system that uses a different two-particle interaction $\hat{W}^{\prime}$, under the mild restriction that the density has to be analytic in timevan Leeuwen (2001). If we choose the two-particle interaction in this other system to be $\hat{W}^{\prime}=0$, the theorem guarantees the existence of the effective potential $v_{s}({\mathbf{r}},t)$ for a Kohn-Sham system that reproduces the same time-dependent density as the interacting system of interest. Here we point out the modifications to the original theorem in order to make it compatible with imaginary-time evolution.

A complex $t$ value does not pose any problems with the operations performed in the original proof, where $t$ appears in some time derivatives but otherwise is treated as a parameter. We add time-dependent uniform potentials $\lambda(t)$ and $\lambda^{\prime}(t)$ to the unprimed and primed Hamiltonians to conserve the norm of the wavefunctions. The origin of these terms will be discussed in the next section. The Hamiltonian $\hat{H}$ of a finite many-particle system is then given by

[TABLE]

expressed in terms of creation and annihilation operators

[TABLE]

Since $\lambda(t)$ is not an operator, it does not affect any of the commutators involving $\hat{H}(t)$ in the various Heisenberg equations of motion underpinning the proof of the van Leeuwen theorem. There is only one detail to note, regarding the freedom to add an arbitrary $C(t)$ to the potential of the primed system, $v^{\prime}({\mathbf{r}},t)$, in the original proof. From Eq. (17b) a time-dependent constant in the potential modifies the Hamiltonian by an additional term $C(t)\hat{N}$, where $\hat{N}$ is the number operator. For the systems of interest, the number of particles $N$ is fixed so $C(t)N$ is a time-dependent uniform potential like $\lambda^{\prime}(t)$, which means that a norm-conserving $\lambda^{\prime}(t)$ will cancel any effect from the choice of $C(t)$. Thus, with $\lambda(t)$ and $\lambda^{\prime}(t)$ chosen to ensure that the norm of states in both the unprimed and primed systems is held at unity, the van Leeuwen theorem holds in imaginary time. This is a powerful result since it allows us to think about imaginary-time propagation in the Kohn-Sham system in terms of what it does in the real system, allowing the Wick-rotation connections from quantum mechanics to statistical mechanics to be employed. For example, it justifies the use of the Kohn-Sham system as a stand-in for the interacting system in our calculations performed for imaginary time path integralsKolesov et al. (2018).

III.2 Maintaining Orthonormalization

Orthonormalization of the single-particle states is equivalent to adding a purely time-dependent function $\lambda(t)$ to the many-body Hamiltonian. This takes care of holding the wavefunction normalized, both in the interacting and Kohn-Sham systems, as well as accounting for the orthogonalization step we use in the Kohn-Sham state propagation.

We first consider the interacting system. In real time propagation, the choice of $\lambda(t)$ does not affect the dynamics of density since this spatially-constant offset in energy only results in changing the phase of the wavefunction:

[TABLE]

In imaginary-time propagation, $\lambda(\tau)$ modifies the imaginary-time propagator $\hat{\mathcal{U}}\quantity(\tau,\tau_{0})$ by a time dependent magnitude,

[TABLE]

If $\lambda(\tau)$ is arbitrary, the norm of the wavefunction will change in time, incorrectly scaling the expectation values of observables like density and energy. The norm of the wavefunction can be held fixed by choosing $\lambda(\tau)$ to counteract the norm-altering effect of $\hat{\mathcal{T}}_{\tau}\exp\quantity(-\int_{\tau_{0}}^{\tau}\hat{H}(\tau^{\prime})\differential{\tau^{\prime}})$ when it acts on $\ket{\Psi(\tau_{0})}$. Note that such a $\lambda(\tau)$ will also depend on the starting state. For example, in the time-independent Hamiltonian case presented in Section II.1, from Eq. (4)

[TABLE]

which implies

[TABLE]

that is, to keep the wavefunction normalized, $\lambda(\tau)$ is such that the energies of the Hamiltonian are measured relative to $\expectationvalue{E(\tau)}$. This result holds more generally for time-dependent Hamiltonians as well, which can be shown by using $\hat{\mathcal{U}}\quantity(\tau,\tau_{0})$ from Eq. (19a) and differentiating the norm-conserving equation $1=\bra{\Psi(\tau_{0})}\hat{\mathcal{U}}^{{\dagger}}\quantity(\tau,\tau_{0})\hat{\mathcal{U}}\quantity(\tau,\tau_{0})\ket{\Psi(\tau_{0})}$ to solve for $\lambda(\tau)$. We will assume that such a $\lambda(\tau)$ is used in the interacting system so that the system always remains normalized.

In the Kohn-Sham system the propagator is given by

[TABLE]

where $\hat{H}_{\textrm{KS}}$ acts on the entire Kohn-Sham many-body wavefunction $\ket{\Phi}$ through its constituent single-particle states $\ket{\phi_{j}}$, see Eq. (5). In general $\lambda_{\textrm{KS}}\quantity(\tau)$ differs from the constant $\lambda(\tau)$ of the interacting system, and in addition to normalizing the many-body state, it can account for orthonormalization of the constituent single-particle states.

Orthonormalization of the occupied single-particle states is an invertible transformation as it preserves the subspace spanned by these linearly-independent states. Representing the orthonormalization by matrix $\mathbf{S}$ and given a single-particle Slater determinant wavefunction $\Phi\quantity(\mathbf{r}_{1},\mathbf{r}_{2},\ldots,\mathbf{r}_{N})$, it is straightforward to show that orthonormalization results in $\tilde{\Phi}\quantity({\mathbf{r}}_{1},{\mathbf{r}}_{2},\ldots,{\mathbf{r}}_{N})=\Phi\quantity({\mathbf{r}}_{1},{\mathbf{r}}_{2},\ldots,{\mathbf{r}}_{N})\det\mathbf{S}$. Thus, the orthonormalization step merely amounts to changing the phase and rescaling the many-body wavefunction.

At the starting time $\tau_{0}$, we assume the Kohn-Sham wavefunction is properly normalized. Following the application of the imaginary-time propagator up to a particular time $\tau$, we represent a particular orthonormalization of the single-particle states by an invertible transformation ${\mathbf{S}}(\tau)$. In order for $\lambda_{\textrm{KS}}(\tau)$ to act like the orthonormalization procedure, we require that

[TABLE]

Note that $\absolutevalue{\det S(\tau)}$ will be continuous since it is the reciprocal of the norm of the unnormalized propagated wavefunction. The phase of $\det S(\tau)$ is not important since it changes the phase of the wavefunction, which will not affect the density. We can therefore use any orthonormalization procedure at each time-step without concern about the continuity of the phase, and a purely real $\lambda_{\textrm{KS}}(\tau)$ satisfying $\hat{\mathcal{U}}_{\lambda_{\textrm{KS}}}(\tau,\tau_{0})=\absolutevalue{\det\mathbf{S}(\tau)}$ for all $\tau>\tau_{0}$ is guaranteed to exist.

The above definitions for norm-conserving $\lambda(\tau)$ and $\lambda_{\textrm{KS}}(\tau)$ conclude the proof of the imaginary time extension to the van Leeuwen theorem presented in Section III.1.

III.3 Monotonically Decreasing Energy

In the Kohn-Sham system the Hamiltonian depends on the density, and thus will in general have eigenenergies and eigenvalues that depend on time. In particular, for the density at time $\tau_{\ell}$, $n({\mathbf{r}},\tau_{\ell})$, we are considering a quantum system with the Hamiltonian $\hat{H}_{\textrm{KS}}[n({\mathbf{r}},\tau_{\ell})]$. By propagating the state of interest in imaginary time using this instantaneous Hamiltonian, we are amplifying the low-energy eigenstates of the current Hamiltonian $\hat{H}_{\textrm{KS}}[n(\tau_{\ell})]$, which in general are different than the low-energy eigenstates of the new Hamiltonian, $\hat{H}_{\textrm{KS}}[n(\tau_{\ell+1})]$, and the resultant state could have a higher energy than the previous state. A good example of this is the commonly-observed divergence of SCF loops without a mixing scheme: the $N$-lowest eigenstates of the Hamiltonian $\hat{H}_{\textrm{KS}}[n_{i}]$ are directly used to compute the next density $n_{i+1}$. This also reveals an interesting limiting case of it-TDDFT. If a KS state is propagated to infinite imaginary time before the density used in the instantaneous Hamiltonian $\hat{H}_{\textrm{KS}}[n(\tau_{\ell})]$ is updated, the propagated state will become the ground state of the present Hamiltonian, which is equivalent to populating the $N$-lowest eigenstates of $\hat{H}_{\textrm{KS}}[n(\tau_{\ell})]$. In this way basic SCF can be thought of as it-TDDFT with infinitely large time-steps when using explicit propagation. Indeed, if the time-step in it-TDDFT is taken to be too large, the total energy will diverge, just like in SCF performed without a mixing scheme.

With a reasonable time-step (usually around $1\hbar/\textrm{Ry}$ or smaller), it-TDDFT monotonically decreases the total energy of the system. The van Leeuwen theorem, which connects the KS system to the interacting system, provides the theoretical backbone for this result. While propagation of the Kohn-Sham system is complicated by the dependence on density, in the true interacting system the evolution in imaginary time has the simple form given in Eq. (4).

III.4 it-TDDFT as an Alternative Theoretical Foundation for Stationary States in DFT

The first step in the majority of DFT calculations is to find a density corresponding to a stationary state. A stationary state is an eigenstate of the Hamiltonian, or equivalently, a state that only changes by a phase when evolved in real time or by a multiplicative factor when evolved in imaginary time. Only the first definition is used in KS systems, as it is implicitly assumed by SCF schemes. In systems that are difficult to converge with SCF, owing to their size or metallic character, the second definition becomes more useful, and it can be applied through the it-TDDFT method.

The KS equations are set as an eigenvalue problem, and thus use the first definition. Once a density $n({\mathbf{r}})$ is found such that a choice of $N$ of the single-particle eigenfunctions $\phi_{j}({\mathbf{r}})$ reproduces the same density through Eq. (9), a stationary state has been determined. SCF is used to find ground states, where the $N$ lowest-energy eigenstates are chosen, and $\Delta$SCFKowalczyk et al. (2011), where a different selection of $N$ single-particle eigenstates is chosen, is used to find excited states. For small systems, insulating systems, and systems with low degeneracy of single-particle states, after a few steps of SCF, the eigenstates rarely change order when sorted by energy from one step to the next. This means that occupied single-particle states have similar character to those from the step before, so the density does not change drastically. In these cases SCF converges well so using the eigenstate definition of stationary states is sensible. However, in large systems and in metallic systems, or if an excited state is desired, the above conditions might not hold, leading to charge sloshing. In principle, the KS equations can still be used to verify a stationary state if the density is perfectly converged. In practice, this definition is inadequate in these difficult systems since a suitable approximate density could appear to be far from convergence if the wrong KS eigenstates are occupied, due to the next SCF step returning a very different density from the one given. In addition, this makes it challenging to determine the quality of a non-converging density. For example, in Section IV we examine the performance of SCF on a ruthenium nanocluster, where we show that some non-converged densities give a reasonable energy estimate for the ground state, while others are incorrect.

To address convergence issues, DFT calculations of metallic systems and systems with high single-particle energy degeneracy are often performed with electronic smearingRabuck and Scuseria (1999), where states near the Fermi level are given fractional occupations to simulate nonzero electronic temperature. This mitigates the problem by ensuring that states near each other in energy have similar fractional occupation. Smearing adds an entropic contribution to the energy, so a balance between obtaining an accurate energy and ease of convergence has to be struck. Electronic smearing is a computational tool and not intended to be an accurate representation of the effects of temperature, so it should be incrementally reduced until the solution with no smearing is achievedMichelini et al. (1998). In fact, cases have been found where even small amounts of electronic smearing produce significantly different results from the same calculation performed with integer occupations, such as a HOMO-LUMO gap energy that differ by one order of magnitudeBasiuk (2011). As we show in the ruthenium nanocluster system in Section IV, achieving convergence while applying electronic smearing can still require finesse and guesswork.

In systems where SCF convergence is hard to attain, instead of using the KS equations to define a stationary state, we can use a state’s invariance under imaginary-time evolution. If a KS wavefunction stays constant when propagated in imaginary time, then its single-particle states span the same subspace as a set of eigenstates which solve the KS equations. The converse is true as well, namely that a set of $N$ KS eigenstates satisfying the KS equations self-consistently will be invariant under imaginary-time evolution, ignoring the possibly changing norm which can be corrected (as discussed in Sec. III.2). Thus, finding a KS many-body state $\ket{\Phi(\tau_{0})}$ such that $\ket{\Phi(\tau)}=\hat{\mathcal{U}}(\tau,\tau_{0})\ket{\Phi\quantity(\tau_{0})}=\ket{\Phi(\tau_{0})}$, where the Hamiltonian in the propagator $\hat{\mathcal{U}}$ contains orthonormality-preserving $\lambda(\tau)$, is equivalent to finding a set of $N$ single-particle states that satisfy the KS equations.

This definition has a few advantages. In systems where the single-particle states are close in energy, occupation ambiguities and charge-sloshing issues are eliminated because it-TDDFT follows the occupied orbitals throughout their evolution. Additionally, it-TDDFT handles systems with degenerate states well since an initial state will converge to one of the states within the degenerate stationary-state subspace without being affected by the unoccupied states of identical energy.

III.5 Practical Advantages of it-TDDFT

One convenience afforded by it-TDDFT is that a user only needs to choose a single parameter, the time-step, when attempting to converge a system. Compare this to the various parameters usually required for SCF with a mixing scheme: the number of past states to mix, the mixing weight, and the amount of electronic smearing, to name a few. When encountering a set of nonconvergent parameters, it is often unclear which direction to change each parameter for a better chance at convergence. In addition, there are systems where different stationary states can be obtained for slight variations in the mixing parameters, as shown in the case of a Cu13 cluster in Section IV. In contrast, convergence in the it-TDDFT method is not very sensitive to the choice of time-step, and any choice smaller than a convergent time-step will lead to the same density trajectory in imaginary time given the same starting state. This property allows us to eliminate this parameter choice if desired through algorithms that automatically adjust the time-step on the fly. We found that the simple procedure of increasing the time-step while total energy decreases, and decreasing the time-step when it does not, can perform nearly as well as using a static convergent time-step that is as large as possible.

Another practical advantage of using imaginary-time evolution is that not-yet-converged states still have physical meaning. The single-particle states and the electronic density used in the KS Hamiltonian are self-consistent at all times, and in principle this density trajectory is equal to the imaginary-time evolving density of the interacting system by the van Leeuwen theorem. Through this connection, the partially-converged KS state corresponds to a superposition of a dominant ground state component and a few low-amplitude excited states. As such, even before the it-TDDFT ground state calculation has converged according to user-specified energy or density tolerance criteria, approximate ground state observables can be computed. This property allows for preliminary calculations of band structure, energies, optical properties, or atomic forces while the ground state calculation continues to be refined. In contrast, there are no guarantees of validity for observables calculated from intermediate states produced in a SCF loop since they are not self-consistent and can be far from the correct KS ground state.

When an SCF loop ceases to make progress, not much is gained aside from the knowledge that the particular set of mixing parameters did not lead to convergence. It could take a subtantial amount of tweaking of these parameters before chancing upon a set that works, consuming time and computational resources. This highlights another strength of it-TDDFT: the calculation time used will always improve the quality of the state at hand. It is also straightforward to continue a calculation from the last saved state, enabling incremental improvement of an approximate stationary state over multiple runs.

In our discussion we have assumed that we are performing DFT with a Kohn-Sham system, which uses a single Slater determinant. It is possible to apply it-TDDFT for finding stationary states in other approaches which use linear combinations of Slater determinants, or ensemble DFT, since a DFT model that reproduces the same density trajectory as the true interacting system will evolve an arbitrary starting state into a stationary state when propagated in imaginary time.

IV Example Calculations

In order to compare two different densities, it is useful to have a measure of distance. We will use half the $L^{1}$ distance for its intuitive physical meaning:

[TABLE]

which can be interpreted as the number of electrons in the wrong place relative to the reference density $n_{0}$. This can be seen by using the fact that both densities integrate to the same value, the total number of electrons. The integral of the absolute value of the density difference over all space adds up the excess density and the negative of the lacking density, both contributing equally, so the $1/2$ factor is needed to obtain the number of electrons out of place.

As a demonstration of using it-TDDFT to determine a ground state, we apply the method to a benzene molecule and show that it produces the same density and energy as a standard SCF calculation. We initialize the single-particle electronic states by drawing basis coefficients from a uniform distribution and orthonormalizing the single-particle wavefunctions. Propagating this initial state in imaginary time, we obtain the same ground state as that determined by an SCF approach with Pulay mixing. The Kohn-Sham total energy $E$ and the density distance $D[n,n_{0}]$ of the propagated state are plotted as a function of imaginary time in Fig. 1, both relative to the SCF-determined ground state. These quantities tend to zero, showing that it-TDDFT indeed produces a Kohn-Sham state that has the same energy and density as the ground state determined with an SCF approach. As an additional check, running SCF at the end of the imaginary-time propagation produces SCF convergence after the first step.

As our next example we consider the Cu13 nanocluster. Hoyt et al.Hoyt et al. (2018) simulated this system in its ground state magnetization of $m=5\mu_{B}$ and in an excited state with magnetization $m=3\mu_{B}$, commenting that the $m=1\mu_{B}$ excited state was tricky to converge, making it a good candidate for our it-TDDFT method. In Fig. 2 and 3, we present the main results of our computations for the self-consistent KS states with $m=1\mu_{B}$ magnetization, which has total spin $1/2$. Additional information can be found in Table 2. SCF has trouble finding the minimum energy states in this fixed-spin system, due to the fact that there are five degenerate statesHoyt et al. (2018). Fig. 2 shows the energy trajectories in imaginary time of 12 different random starting configurations, and each of these converges to one of five lowest-energy states. To help identify the equality of final states, for both the it-TDDFT and SCF runs, we also computed the density distances between each combination of obtained states. States with energies within $10^{-2}\,\textrm{meV}$ and a density distance of less than $(1/100)e$ of each other were considered equal. The ground state of the system, with magnetization $m=5\mu_{B}$, contains five degenerate valence electrons which have unpaired spinsHoyt et al. (2018). To obtain a magnetization of $m=1\mu_{B}$, four of the electrons need to pair spins, leaving five possibilities of which electron remains unpaired. Fig. 2 shows the magnetization density, defined as the difference between spin up and spin down electron density, of one of these five lowest-energy states. There are five equivalent ways to place such a magnetization density on the icosahedral shape of the copper cluster, explaining the degeneracy. The visible differences of up to $0.1\,\textrm{meV}$ in the energies of these degenerate states are due to the discretization effects of the real space grid breaking the icosahedral symmetry.

Our approach is better at finding the lowest-energy states compared to SCF, which for different mixing parameters often converges to other excited states of spin $+1/2$ (the energies of which are shown as dashed lines in the figure). In Fig. 3, we show the electronic energies of the states obtained using SCF with Pulay mixing, for various mixing parameter choices. Even small changes in the mixing parameters can result in a different final state. This happens in metallic systems where the gap between occupied and unoccupied states is small, causing SCF to find a low-lying excited state.

For our final example of applying our method we consider the Ru55 nanocluster. Montemore et al.Montemore et al. (2018) studied catalysis on the surface of this structure, and found that the spin-unpolarized ground state calculation was difficult to converge with SCF.

In Table 1, we show the results of using SCF with Pulay mixing to find the ground state of the spin-unpolarized Ru55 cluster. The number of past densities to mix was kept at $n=5$ for all trials and mixing weights ranging from $0.02$ to $0.20$ were tested. We used Fermi electronic smearing for half the trials with $T=300\,\textrm{K}$. For each run, we list the energy of the final step relative to the energy calculated with it-TDDFT, the density difference $\Delta\rho_{\textrm{max}}$, and whether the run converged or not. The density difference $\Delta\rho_{\textrm{max}}$ refers to the maximum elementwise difference in the density matrix between the final and penultimate step and is typically used to determine convergence. We used the criterion $\Delta\rho_{\textrm{max}}<10^{-6}$. Only a few mixing weights result in convergence, namely the smallest ones with $T=300\,\textrm{K}$ of smearing. In these runs, the entropic energy contribution is $78\,\textrm{meV}$ relative to the ground state energy. None of the runs without electronic smearing converge, despite the fact that some parameter configurations obtain energies similar to the ground state energy. The states resulting from unconverged runs generally should not be trusted as they may not be acceptable approximations to actual solutions, which have to satisfy the KS equations self-consistently. For example, in Table 1, examining the row with mixing weight $0.10$ and comparing the $T=0\,\textrm{K}$ and $T=300\,\textrm{K}$ cases, we find that even though $\Delta\rho_{\textrm{max}}=0.589$ in the latter run is smaller than $\Delta\rho_{\textrm{max}}=0.702$ in the former, the energy of the state obtained with $T=300\,\textrm{K}$ is more than $6\,\textrm{eV}$ off from the correct ground state energy while the $T=0\,\textrm{K}$ run is only about $0.006\,\textrm{eV}$ off. Applying our it-TDDFT method to the Ru55 cluster produces the ground state without issue, as illustrated in Fig. 4. The observed monotonically decreasing energy and density distance $D[n,n_{0}]$ show consistent progress, as we expect from the theory.

V Conclusion

The first step of any Kohn-Sham DFT calculation is the determination of a self-consistent solution to the KS equations, resulting in a density corresponding to a stationary state of the many-body interacting system. While the standard method of using the iterative SCF procedure generally produces a solution efficiently, there are important classes of systems that pose problems for this approach due to their small band gaps or degenerate single-particle energies. We have proposed the it-TDDFT method as an alternative means for solving the KS equations in these difficult systems, and shown how it avoids the issues which affect SCF.

We established that the van Leeuwen theorem, a key theoretical foundation for TDDFT methods, can be extended to imaginary time, thereby ensuring convergence to a stationary state independent of the exchange-correlation potential and level of theory used in the model system. In addition, we discussed how it-TDDFT could be used in an alternative but equivalent definition of stationary states in DFT, better suited for metallic systems and systems with degenerate or nearly-degenerate states and based on the time-dependent Kohn-Sham equations. The it-TDDFT method also exhibits a number of practical advantages, such as justifying approximations to observables of interest before the ground state calculation is fully converged, requiring few input parameters, and allowing easy refinements of the results of previous runs by continuing from a saved state.

In the copper and ruthenium nanoclusters considered here, we demonstrated how SCF can struggle to find the electronic ground state, either converging to low-lying excited states or getting stuck in charge-sloshing cycles. These systems were readily converged by it-TDDFT, showcasing its robustness through smooth trajectories with monotonically decreasing energy. For these systems we either ran the calculation as spin-unpolarized, or with a fixed total spin. This is not an inherent limitation of the method, as one could simply run the calculation with all possible spin polarizations and select the state with the lowest energy. The method can be adapted to non-collinear spin systems, since the operating principle depends only on the Hamiltonian being able to differentiate states by energy. Furthermore, while we used finite systems for our example calculations, our method can be extended to find ground states of periodic systems by simultaneously propagating Kohn-Sham states at multiple $k$-points.

Given an existing TDDFT code which evolves systems in real-time, it should be relatively straightforward to implement a prototype of the presented it-TDDFT approach, requiring only an imaginary time substitution in the propagation step and a method to orthonormalize the single-particle states. While more efficient implementations could be examined in the future, the low barrier to utilizing it-TDDFT could make it an attractive alternative option for those dealing with particularly vexing systems.

VI Supplemental Materials