Accurate Evaluation of Real-Time Density Functional Theory Providing

Subscriber access provided by UNIV OF DURHAM

Article

A Different Approach to the Accurate Evaluation of Real-Time Density Functional Theory Providing Access to Challenging Electron Dynamics Ingo Schelter, and Stephan Kummel J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01013 • Publication Date (Web): 28 Feb 2018 Downloaded from http://pubs.acs.org on March 4, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

A Different Approach to the Accurate Evaluation of Real-Time Density Functional Theory Providing Access to Challenging Electron Dynamics Ingo Schelter∗ and Stephan Kümmel∗ University of Bayreuth, Germany E-mail: [email protected]; [email protected]

Abstract We demonstrate that electronic excitations and their transition densities can be obtained from real-time density functional theory calculations with great accuracy by relating the data from numerical propagation to the analytical form of the electronic response after a boost excitation. The latter is derived in this article. This approach facilitates to obtain oscillator strengths quantitatively, to identify excitations that carry very small oscillator strengths, and to assess electronic couplings from transition densities based on comparatively short propagation times. These features are of interest in particular for studying light-harvesting systems. For demonstration purposes we study the Q-band excitations of bacteriochlorophyll a (BChla) and calculate coupling strengths between two BChla to check the validity of the dipole-dipole and pure Coulomb coupling mechanisms. For further illustration we investigate the paradigm test systems Na4 and the coupling between two Na2 dimers.

1

Introduction

These aspects come together in a particularly challenging way in the research on natural light harvesting systems. The photosynthetic process transforms solar energy delivered in the form of photons into other forms of energy with such a robustness and efficiency that its microscopic understanding, and its mimicking in man-made devices, are considered to be one of today’s greatest and most relevant intellectual challenges. The highly efficient absorption of light and the subsequent transfer of excitation energy in the complex molecular aggregates that form the photosynthetic complexes of plants and bacteria constitutes a very formidable quantum mechanical many-particle problem. Due to the large number of particles – typically thousands of nuclei and electrons – that are involved, a first principles description based on correlated

Density functional theory (DFT) 1 and time-dependent density functional theory 2–4 (TDDFT) have emerged to be among the most frequently used many-particle theories for calculating electronic structure properties and electron dynamics. To a large extent this is a consequence of the favorable ratio of accuracy to computational cost which can be achieved with (TD)DFT in many situations of practical interest. 5–14 This practical usefulness, however, hinges on the one hand on the availability of functionals that reliably approximate the many-particle effects of exchange and correlation (xc), and on the other hand on computational methods which realize such functionals with an efficiency that is high enough to allow studies on systems of practical relevance.

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

wave functions will be impossible for years to come. However, also (TD)DFT is seriously challenged by the size of the problem, and methods that scale well with respect to system size must be brought to bear. The real-time approach to TDDFT, 15,16 which gives access to linear and non-linear phenomena, is ideally suited for accepting this challenge. It scales favorably with the system size because adding an electron just adds one more orbital to propagate and no unoccupied states are needed. The corresponding equations are also easy to parallelize as the time step for each orbital hardly needs information about other orbitals, i.e., communication needs are intrinsically small. On the technical side, there have been two development strands. On the one hand, real-time TDDFT started with implementations that represent the orbitals on real-space grids. 5,7–9,15,17–19 In some cases, the decision for a grid has been related to the fact that grids represent orbitals without a particular focus on specific regions of space, contrary to atom-centered basis sets, and are therefore well suited for describing non-perturbative dynamics 7,18,20–23 and long-range emission processes. 24–26 However, real-space grids have also remained popular because of the excellent further parallelization that they offer due to their underlying sparse matrix algebra, leading to efficient codes with many applications. 11,27–31 On the other hand, real-time TDDFT is also increasingly used in combination with basis sets. 14,16,32–40 In particular atom centered basis set codes can profit from great efficiency when bases of moderate size can be used. Furthermore, they are particularly helpful for functionals that require the evaluation of Fock integrals such as hybrids, 16 as the latter are expensive to compute on real-space grids. A possible drawback of the propagation technique may be seen in the fact that such calculations return, at least directly, just the time-dependent density n(r, t), whereas the linearized Kohn-Sham equations, 2,41–43 as used, e.g., in the Casida formalism, directly return excitation energies and oscillator strengths. I.e., they return information on dipole forbidden

Page 2 of 30

excitations or ones with very small oscillator strengths. Furthermore, it has been noted early 18 that obtaining spectral observables from real-time propagation needs extensive propagation times for a satisfying spectral resolution, unless special measures are taken as discussed in detail below. Therefore, it can be quite hard to identify excitations carrying little oscillator strength as well as excitations that lie energetically close to nearby stronger ones in the realtime approach. These challenges become particularly relevant when studying organic lightharvesting systems, as these typically contain very many electrons so that long propagation times come at an immense computational cost, yet they also show charge-transfer (CT) and thus almost dark excitations as well as closelying excitations, making a high spectral resolution mandatory. In recent years, techniques have been developed which aim at extracting the electric response from real-time simulations with increased accuracy and reduced computational cost. For example, in analogy to techniques used in static response calculations, 44,45 fitting procedures have been developed to extract response coefficients. 30,38 However, they require several propagation runs with perturbations of different frequencies and/or field strengths. Thus, they are still computationally demanding and better suited for obtaining non-linear response coefficients of small molecules than for obtaining insight into the reponse of large systems. Different approaches for increasing computational efficiency have been put forward. Andrade et al. 11 suggested a compressed sensing strategy. Alternatively, Bruner et al. 14 showed that the propagation time required for obtaining the photoabsorption spectrum can be considerably reduced by fitting the dipole signal from a relatively short propagation with a Padé approximant. In this article we take a look at the same problem – extracting the response accurately from short propagation runs – from a different perspective. By fitting the analytical form of the dipole signal after a boost excitation to the numerical dipole signal in frequency space, we obtain excitation energies and oscillator strengths


2

Page 3 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with great accuracy. We further show how transition densities can be obtained accurately from relatively short propagation times. Our techniques are very helpful when studying large systems. Our first tests are done for the small Na4 cluster, which is an established TDDFT benchmark system, and for a bacteriochlorophyll a (BChla) molecule, which is a common chromophore in bacterial light-harvesting systems. As the coupling strength between chromophores is an observable of key interest in questions of energy transfer, we also evaluate the coupling between pairs of these molecules. In doing so we compare Förster-type dipole-dipole coupling 46 to pure Coulomb coupling between the full transition densities 46,47 and to the coupling strength that one infers from the Davydov splitting. 48,49 We find that the error of the dipoledipole approximation falls below the 10%-limit for separations larger than 20 Bohr radii (a0 ) in the case of the two Na2 molecules, whereas two BChla must be separated by at least 50 a0 for the dipole-dipole coupling to become a reasonable (10%-error) approximation. The article is organized as follows. First, we briefly review real-time Kohn-Sham (KS) TDDFT and the simplest evaluation of electronic spectra and transition densities in Section 2. In Sections 3 and 4 we explain our new method for evaluating spectra and transition densities and demonstrate that we can thus obtain excitations of Na4 and BChla that remain obscure when the real-time density is evaluated in the simple way. The calculation of coupling strengths is discussed in Section 5 and we offer conclusions in Section 6. Numerical details, in particular about the Bayreuth-DFT (BTDFT) program that we developed and used for this work, are given in appendices.

2

our work and our notation, we briefly reiterate the basic equations.

2.1

Time-Dependent Kohn-Sham Scheme

The KS formulation 52 is the most common form of TDDFT. The time-dependent KS equations i~

∂ ˆ KS (r, t) ϕj (r, t) ϕj (r, t) = h ∂t

(1)

describe the time evolution of the KS orbitals ϕj . In eq (1) 2 ˆ KS (r, t) = − ~ ∇2 + vKS (r, t) h 2m

(2)

is the KS Hamiltonian, where ~ is Plank’s constant, m is the electron mass, and i the imaginary unit. The spin index is suppressed for simplicity. The local KS potential vKS (r, t) = vext (r, t) + vH (r, t) + vxc (r, t) (3) consists of an external potential vext that describes external fields and the nuclear potentials, the Hartree potential vH , and the xc potential vxc . vH and vxc are functionals of the electron density which is expressed by the KS orbitals n (r, t) =

N occ X j=1

|ϕj (r, t)|2

(4)

with Nocc being the number of occupied orbitals. The approximation used for the unknown exact vxc is decisive for the accuracy of a (TD)DFT calculation. The adiabatic local density approximation 51 (TDLDA) is widely used in real-time TDDFT. 7,8,53 In addition to a low computational cost it offers the advantage of being robust, i.e., its shortcomings are well known and predictable. We use it here (in the parametrization of Ref. 54) because our aim is to demonstrate the efficient evaluation of realtime observables for transparent test cases. It is well known, though, that TDLDA has systematic shortcomings, e.g., with respect to CT

Basic Equations of TimeDependent Density Functional Theory in RealTime

TDDFT is well known and established since many years. 2–4,16,50,51 To clarify the context of


3


excitations, that can be remedied with functionals of greater complexity. 55 An investigation of whether or which properties of lightharvesting systems can be computed reasonably with TDLDA is the topic of a different study. 56 Since eq (1) is a differential equation in space and time, boundary and initial conditions are required for its solution. For finite systems, zero boundary conditions typically apply, i.e., the orbitals approach zero at infinity (in our case on the boundary of a numerical grid). An exception to this rule are ionization processes in which density escapes “to infinity”, 57,58 but such processes are not of interest in the present context. The initial state has to be compatible with the Runge-Gross theorem 51 and is typically the ground state.

2.2

step t → t + ∆t is taken in the approximation i

= τˆe

Rt

t0

ˆ KS dt′ h

(5)

ϕj (r, t0 )

(7)

where the system’s state is known at time t and ˆ KS is taken at t+ ∆t . ∆t has to be small enough h 2 for the KS Hamiltonian to be approximately constant during the time step. Since the KS Hamiltonian depends on the time-dependent density through the KS potential vKS , one ˆ KS (t + ∆t ) by a predictor-corrector obtains h 2 scheme and/or by extrapolating vKS . 19,39,62 For consistency and stability reasons we always used a predictor-corrector scheme for the calculations presented in this article. i ˆ By applying e 2~ hKS ∆t to eq (7) from the left and expanding the exponential into its Taylor series up to first order in ∆t one arrives at the Crank-Nicolson propagator 62 i ˆ 1 + hKS ∆t ϕi (r, t + ∆t) = 2~ (8) i ˆ 1 − hKS ∆t ϕi (r, t) . 2~

In order to solve eq (1) with the real-time approach it is practical to use the propagator Uˆ (t, t0 )

− ~i

ˆ

ϕi (r, t + ∆t) = e− ~ hKS ∆t ϕi (r, t) ,

Real-Time Propagation

ϕj (r, t) = Uˆ (t, t0 )ϕj (r, t0 )

Page 4 of 30

The Crank-Nicolson propagator requires the solution of a differential equation in r for ϕj (r, t+ ∆t). The latter turns into a sparse system of linear equations on the equidistant real-space grid that we use. The Crank-Nicolson propagator is norm conserving and unconditionally stable for a given KS potential. 62 The latter means that there is no additional constraint on ∆t other than the one of being small compared to the time scales of the observed dynamics of the system. However, in practice we found, similar to Repisky et al. 39 for the second-order Magnus propagator, that a propagation without a predictor-corrector scheme is prone to instabilities. In our parallelized grid-based code we compared the Crank-Nicolson propagator to a fourth-order Taylor propagator. 19,60 We found that typically the Crank-Nicolson propagator is more efficient, because it allows for a larger time step and the resulting smaller total number of time steps usually overcompensates its higher computation cost per step.

(6)

and discretize the time coordinate by using a finite time step ∆t. In typical propagation schemes, the action of the time ordering operator τˆ is included per construction. There are numerous ways of how the propagator can be represented numerically. As the intention of our article is not a review or discussion of propagator approximations, we refer the reader to Refs. 7,19,30,31,33,35–37,40,59–61 for different ways of how efficient numerical propagation can be achieved. We stress that the method that our article is about, i.e., the efficient evaluation of spectral observables, can be combined with any propagator as it just uses the timedependent density. In the following, we therefore only briefly summarize how our own calculations are technically done. In our calculations, quantities such as orbitals and potentials are represented on a real-space grid. In the notation of the following equations, we omit the discretization of the spatial coordinates for the sake of clarity. A propagation


4

Page 5 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.3

Elementary Evaluation Electronic Excitations

ened and consequently the line shapes in the frequency domain are not δ-function-like peaks, but oscillate spuriously. Sometimes the dipole moment is damped artificially, e.g., by a decaying exponential function, such that one obtains a smooth spectrum without the spurious oscillations. This procedure, however, is accompanied by a loss of spectral resolution. 5,18,39 A second quantity that is often used to identify excitations is the PSD 18,60,62,63

of

In real-time TDDFT electronic excitations are traditionally identified from a dipole strength function 4,5,17 (DSF) or a power spectral density 18,60,62 (PSD), which are extracted from the time-dependent dipole moment after an excitation. In this subsection we recapitulate these evaluation techniques. Following Yabana and Bertsch 5,17 the electronic ground state of the finite system with 0 ground-state orbitals ϕj is excited by a boost at time t0 = 0 to obtain the initial state for the propagation, i.e., ϕj (r, t0 = 0) = eir·k ϕ0j (r) .

P (~ω) =

(ϑ)

(11)

=P (ϑ)

The PSD shows peaks at the same frequencies as the DSF. It has been used as an alternative to the DSF because for some systems the more pronounced peak structure of the PSD allows for a better identification of the excitation energies. The peak heights of the PSD, however, do not directly correspond to the oscillator strengths. In the following we exemplify these evaluation techniques. The first system that we study is the small cluster Na4 , which we chose because sodium clusters are frequently used paradigm test systems. 24,64,65 The PSD for the Na4 cluster resulting from a 50 fs propagation is shown in Figure 1 (a). This kind of representation of the PSD is frequently found in the literature. 3–5,7,17,18,60,66–70 The cluster is oriented such that its atoms lie in the yz-plane with its long axis along z and the short axis along y (cf. inset in Figure 1). From the spectrum up to five spikes indicating electronic excitations at about 1.8 eV, 2.6 eV, 3.0 eV, 3.3 eV, and 3.5 eV can be identified. Their respective energies can be deduced from the peak positions, but only with a relatively poor energy resolution of ≈ 0.1 eV. 60 As a second example we show the PSD of the Q-band 71–73 of BChla from a 100 fs propagation in Figure 2 (a), with the inset depicting the molecular structure (coordinates taken from Oviedo; 74 cf. Appendix F for further details). Two spikes can be seen that correspond to the Qy (at 1.7 eV) and Qx (at 2.0 eV) transitions.

(9)

Here, k determines the direction and energy 2 |k|2 Eboost = N ~ 2m of the boost, where N is the number of electrons. A weak boost corresponds to an instantaneous, dipole-like kick that excites the many-electron system with all frequencies (cf. Section 3). The DSF, which corresponds to the photoabsorption spectrum, is usually extracted as the imaginary part of the Fourier transform of the induced dipole moment δµ(t) = µ(t) − µ0 via 4,17 S(~ω) = −

X 2 1 X δ µ . 2 ˜(ϑ) γ (ω) 3 ϑ=x,y,z γ=x,y,z {z } |

(ϑ) 2mω Im Tr δ µ ˜γ (ω) . (10) 2 3πe~ k

µγ denotes the γ component of the dipole moment after a boost in the ϑ direction with γ, ϑ ∈ {x, y, z}. 4 Ground-state quantities such as the static dipole moment µ0 are in the following indicated by the subscript 0. The tilde marks a Fourier transform f˜(ω) = F [f (t)] and Im denotes taking theh imaginary part. The trace is i (ϑ) (x) (y) (z) evaluated as Tr fγ = fx + fy + fz . After the boost excitation, the dipole moment oscillates with the frequencies that correspond to the energies of the dipole active excitations of the electronic system (cf. Section 3) whose transition dipoles are not orthogonal to the boost direction. The DSF shows peaks at these energies with heights proportional to the corresponding oscillator strengths. Due to the finite propagation time, the peaks are broad-


5


4

a

z

1.2 Power spectrum P(-hω) [arbitrary units]

3.5

5.78 a0

3 2.5 2 1.5 1

y

0.5 0 30 25

Dipole strength S(-hω) [1/eV]

PSD

Refined DSF Fit

b

15 10 5 0

1

c

Evaluated spectrum Along Z Along Y Y2 Along X

Z1

0.8 X2 0.6 0.4

X3

0.2

Z2

Y1

X1

Z3

Z4

Y3

0 1.5

2

2.5 3 Excitation energy [eV]

PSD

b

Refined DSF Fit

0.6 0.4 0.2 0

10 8 6 4 2 0 -2

Oscillator strength (Peak height)

1.2

a

0.8

20

-5 Oscillator strength (Peak height)

1

12 Dipole strength S(-hω) [1/eV]

Power spectrum P(-hω) [arbitrary units]

4.5

12.28 a0

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 30

3.5

0.25

Evaluated spectrum Qy CT1 CT2 Qx

Qy

c

0.2 0.15 0.1

Qx

0.05

CT1

CT2

0 1.5

Figure 1: TDLDA dipole spectrum of the Na4 cluster after a 50 fs propagation showing (a) the power spectral density, (b) the imaginary part of the Fourier transformed dipole moment scaled with Tω and with increased sampling rate (the arrows indicate excitations) and (c) the purified spectrum broadened with Gaussian func2 tions e−[~(ω−ω0j )/η] with η = 0.05 eV. The excitations are named after their symmetry with respect to the coordinate system (and molecular axis) and their energetic ordering.

1.6

1.7 1.8 1.9 Excitation energy [eV]

2

2.1

Figure 2: TDLDA dipole spectrum of BChla after a 100 fs propagation showing (a) the power spectral density, (b) the imaginary part of the Fourier transformed dipole moment scaled with ω and with increased sampling rate (the arT rows indicate excitations) and (c) the purified spectrum broadened with Gaussian functions 2 e−[~(ω−ω0j )/η] with η = 0.025 eV.


6

Page 7 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.4


Elementary Evaluation Transition Densities

lying close to nearby stronger ones. The problem can be resolved by greatly increasing the propagation time, which in turn increases the spectral resolution. Long propagation times, however, severely increase the computational cost. Repisky et al. 39 have reported a scheme for how to decompose transition densities obtained from the four-component Dirac-KohnSham equation into molecular orbital contributions exploiting properties of their Gaussian basis set. Our aim in the following is a different one: We do not decompose the transition density, but show how to calculate transition densities accurately without increasing the propagation time and only using the time-dependent density.

of

Transition densities 75 characterize excitations and allow for the intuitively accessible visualization of excitation dynamics. 63,69,76 In the context of light-harvesting systems they can be particularly useful for calculating Coulomb-like coupling strengths 46 between chromophores in a donor-acceptor model, as we discuss in Section 5. The transition density for the excitation from the many-electron ground state |0i to the excited state |ji with transition energy ~ω0j is defined by ρ0j = hj|ˆ n|0i. 75 In real-time TDDFT it can be obtained from the imaginary part of the Fourier transformed induced density fluctuation δn (r, t) = n(r, t) − n0 (r) evaluated at the desired excitation energy, 69,70,76 ρ0j (r) ∝ Im {δ˜ n (r, ω0j )}

3

(12)

(see further below for a derivation of this relation). From this perspective it becomes clear that the transition density can be interpreted as showing a real-space “snapshot” of how the density oscillates at the excitation energy ~ω0j . While there is per se nothing wrong with calculating a transition density in this way from real-time TDDFT, in practice the approach is often limited in its accuracy and resolution. The reason is related to the above discussed limited energy resolution of the traditional approach. For calculating the transition density one first has to obtain the excitation energy as described above, and then evaluate eq (12) at this energy. Figure 1 (a) and 2 (a) demonstrate that in the traditional approach the excitation energy is obtained only with some uncertainty, which limits the accuracy of the corresponding transition density. Moreover, two transitions lying energetically close to each other overlap spectrally due to the finite line width (which results from the finite propagation time) and can often hardly be resolved. Consequently, their visualized transition densities as evaluated through eq (12) will show a mixture of the two close lying transitions. This is especially critical for weak transitions

Evaluation of Electronic Excitations from Short Propagation Times

The basic idea of our scheme for improving the accuracy of the evaluation of real-time TDDFT calculations is simple: One can go beyond the PSD and DSF as they are traditionally used by fitting the numerical dipole signals with appropriate analytical functions and thus extract excitation energies and oscillator strengths with high accuracy. The in-principle possibility for such a procedure has already been pointed out, 18 but to the best of our knowledge it has not been realized so far. We here demonstrate that such a procedure, when properly worked out, can increase the accuracy substantially. For setting up this accurate analysis of the dipole signals we first derive analytically which form the DSF takes for a real-time calculation that is initialized by a boost. Then we fit the corresponding function to the numerically calculated signal to obtain excitation energies and oscillator strengths accurately. In Appendix H we compare our approach to the above mentioned one by Bruner et al., 14 who approached the same problem from a different perspective via Padé approximants.


7


3.1

Correspondingly, the products appearing in eq (16) simplify to

Analytical Form of the TimeDependent Dipole Moment after a Boost Excitation

c∗l cj ≈ δl0 δj0 + ik · (δl0 rj0 − δj0 rl0 ) .

For the derivation of the analytical form that the time-dependent dipole moment takes after a boost excitation we start from the quantum mechanical many-body problem. Initially, i.e., for t < 0, the system is in the manyelectron ground state |0i. At t = 0 it is subjected to an instantaneous boost excitation 5,7,14,17,18,28,39,63,66,69,70,76,77 |ψ(t = 0)i = e

ik·ˆ r

(18)

The δl0 δj0 term leads to the static dipole moment µ0 = h0| − eˆ r |0i which is of no interest here because it is not related to excitations. The remaining induced time-dependent dipole moment can be written as δµ(t) ≈ −ie

(13)

|0i .

Page 8 of 30

= −ie

Since the Hamiltonian is again timeindependent for t > 0, the time-dependent wave function |ψ(t)i of the many-electron system can be expanded into eigenfunctions {|ji} of the many-body ground-state Hamiltonian with eigenvalues Ej = ~ωj

= −2e

∞ X

j,l=0 ∞ X j=0 ∞ X j=1

k · (δl0 rj0 − δj0 rl0 ) e−iωlj t rlj

(k · r0j ) r0j e−iω0j t − eiω0j t (k · r0j ) r0j sin(ω0j t) ,

(19)

(16)

where we made use of ωjl = −ωlj and rl0 = r0l ∈ R. Thus, the dipole moment is a superposition of harmonic oscillations with the angular frequencies corresponding to the 0 → j transitions. The amplitude of each harmonic oscillation depends on the boost strength k = |k|, the strength of the transition dipole |−er0j | and the angle ϕ0j between boost vector and transition dipole through the inner product (k · r0j ) = |k| |r0j | cos(ϕ0j ). At this point the problem has been reduced to calculating the induced time-dependent dipole moment. The latter can simply be extracted from a real-time TDDFT calculation Z δµ(t) = −e δn (r, t) r d3 r . (20)

with ωlj = ωj − ωl and rlj = hl| rˆ |ji. The asterisk denotes complex conjugation. When the boost excitation is weak, i.e., in the linear response regime, this expression can be simplified further by expanding the exponential eik·ˆr into a Taylor series and neglecting non-linear terms (cf. Appendix A for a discussion). Thus, the expansion coefficients from eq (15) are approximately

The most accurate way to extract dipole-active excitations is through three calculations, each with a boost k(ϑ) = kˆ eϑ along one of three orthogonal directions, e.g., the Cartesian ones (i.e. ϑ ∈ {x, y, z}). eˆϑ denotes the corresponding unit vector. This leads to three (time-dependent) dipole moment vectors δµ(ϑ) (eq (19)) that can be written as a 3 × 3 matrix (ϑ) with components δµγ , where γ ∈ {x, y, z}. (ϑ) The dependence on the angles ϕ0j can be re-

|ψ(t)i =

∞ X j=0

cj |ji e−iωj t ,

(14)

where the {|ji} can be chosen to be real-valued, and the expansion coefficients cj in eq (14) are constant for t > 0 cj = hj |ψ(t = 0)i = hj| eik·ˆr |0i .

(15)

From the time-dependent wave function results a time-dependent dipole moment µ(t) = hψ(t)| − eˆ r |ψ(t)i ∞ X = −e c∗l cj e−iωlj t rlj j,l=0

cj ≈ hj| 1 + ik · rˆ |0i = δj0 + ik · rj0 .

(17)


8

Page 9 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the ω > 0-part of the spectrum. The propagation time T appears as the inverse line width. Inserting eq (23) into the DSF in eq (10) results in

moved by building the trace X (l) Tr δµ(ϑ) = δµl γ l=x,y,z

= −2ek

∞ X

(r0j )2 sin(ω0j t)

j=1 ∞ X

3e~k =− m

j=1

f0j sin(ω0j t) , ω0j

S(~ω) =

f0j

j=1

(21)

T →∞

−→

where

2mω0j |r0j |2 (22) 3~ denotes the usual P oscillator strengths. 75 In eq (21) we used (kˆ e · r )r = l=x,y,z 2 l 0j 0j,l 2 2 2 k r0j,x + r0j,y + r0j,z = k (r0j ) and r0l ∈ R.

3.3

f0j =

3.2

∞ X

∞ X j=0

ω sin((ω − ω0j )T ) ω0j π~(ω − ω0j ) f0j δ(~ω − ~ω0j ) .

(24)

Accurate Excitations from Real-Time Propagation

Equation (23) or (24), respectively, is the decisive relation for our approach to evaluating spectra accurately from relatively short propagation times. To this end, we proceed as follows. We first solve the time-dependent KS equations for three boost directions ϑ ∈ {x, y, z} and as usually obtain n(ϑ) (r, t), compute the time-dependent dipole moment for each boost ϑ from eq (20), and then the DSF from eq (10). The discrete Fourier transform of the dipole moment that is needed in the last step is usually done by some algorithm, e.g., the fast Fourier transform 62 that requires the number of sampling points to be a power of 2. It is thus a common practice 39,62 to append zeros to the discrete time-dependent data up to T (ft) > T such that the upper constraint is fulfilled. The sampling rate in frequency space is then given by ∆ω = T2π (ft) . It is important to note that the spectral lines, i.e., their shapes and positions, remain untouched by adding the zeros and only the sampling rate in frequency space is influenced. Using the next possible value for T (ft) > T and typical propagation times usually leads to the rather spiky spectra as shown in 1 (a) and 2 (a). However, the sampling rate can be improved arbitrarily by choosing a larger T (ft) . The thus refined DSF is drawn in Figures 1 (b) and 2 (b) (label: “Refined DSF”) for Na4 and BChla, respectively. This little numerical trick now allows for the decisive next step: Since we know that the DSF from a finite propagation must have the form of eq (24) we can fit the parameters f0j and ω0j from (24) to the tightly

Analytical Form of the Dipole Strength Function for a Finite Propagation Time

The trace of the dipole moment tensor in eq (21) oscillates with the excitation frequencies f . The Fourier transω0j and amplitudes ∝ ω0j 0j form of eq (21) is purely imaginary and shows sharp, δ-function like peaks at the excitation frequencies, with strengths that are related to the oscillator strengths. However, in practice the propagation time T is always finite and one therefore cannot take the Fourier transform for t from −∞ to +∞, but only within a finite time interval. As shown in Appendix B, this leads to a non-vanishing real part and peaks that are of sin(ωT )/ω shape. As a consequence, the imaginary part of the Fourier transform of eq (21) – which is still the relevant part for photoabsorption – from a finite real-time propagation takes the form (cf. Appendix B) (ϑ) Im Tr δ µ ˜γ (ω) = (23) ∞ 3e~k X f0j sin((ω − ω0j )T ) − , 2m j=1 ω0j ω − ω0j where we neglected the contributions from negative ω0j . This is valid if ω0j T ≫ 1 for all transitions 0 → j, which ensures that the tails of the lines at negative ω0j do not penetrate into


9


sampled, refined DSF line (label: “Fit” in Figures 1 (b) and 2 (b)). In this way, both the excitation energies and the oscillator strengths can be obtained with a greatly improved accuracy. In practice, weak excitations are best identified after fitting stronger excitations and the fit can straightforwardly be done using a publicly available program such as gnuplot. 78 In order to demonstrate that one can obtain considerably better accuracy from a given propagation run with this new evaluation procedure we return to the Na4 cluster, for which five excitations were found from the PSD in 1 (a). With the new procedure eleven excitations can be identified, as shown in 1 (b) and (c). The fit of eq (24) matches the simulation data perfectly. The energies and oscillator strengths shown in 1 (c) were broadened by Gaussian functions to guide the eye. We clearly see that much more information can be extracted from the same propagation run with this evaluation procedure. The practical relevance of the improved accuracy can further be demonstrated by the BChla example. The spectrum of BChla as newly evaluated in Figures 2 (b) and (c) shows the large Qy and Qx excitations that could also be identified in panel (a) of the figure. The fit allows us to assign precise excitation energies and oscillator strengths: The Qy state appears at 1.677 eV with a strength of 0.256, the Qx state at 1.994 eV with a strength of 0.075(1) . Comparing to Ref. 74 from which the molecular geometry has been taken and to the experiment, 79 we see that these excitations are well described by real-time TDLDA. However, in addition to the large Q-excitations, two more, weak excitations become visible in the fit, as clearly seen in panel (c). The latter can be identified as states with a certain CT character that are spuriously shifted down in energy 80,81 when TDLDA is used. It would have been very hard if not impossible to identify these weak excitations, especially the CT1 that lies close to the Qy tran-

Page 10 of 30

sition, from the traditional evaluation shown in panel (a). A comparison with the aforementioned Padé approximant evaluation is given in Appendix H. We close this subsection with three additional remarks. First, in addition to excitation energies and oscillator strengths, the information about the direction of the transition dipoles is often relevant. Examples are the determination of coupling strengths between chromophores or the evaluation of transition densities as discussed in the next sections. While the absolute value of a transition dipole is given by the corresponding oscillator strength via eq (22), one can extract its direction from a single calculation with a boost in a specific direction ϑ. Of course, this procedure is most accurate if boost and respective transition dipole are nearly parallel to each other. To get the transition dipole directions one makes use of eq (19), which states that δ µ(ω ˜ 0j ) ∝ r0j . One now evaluates the sin(ϑ) gle dipole components µγ with γ ∈ {x, y, z} in the way previously explained in this subsection (ϑ) and receives three spectra Sγ with line heights (ϑ) (f0j )γ for each excitation 0 → j. For a specific excitation the ratio of those line heights directly corresponds to the direction of its dipole moment, i.e.,   (ϑ) (f0j )x   (25) r0j ∝ (f0j )(ϑ) y  . (ϑ) (f0j )z This follows from the derivation in this subsection for each dipole component without building the trace in eq (21). Second, in the scheme as outlined above the dependencies on the angles ϕ0j between boost and respective transition dipole from eq (19) are eliminated by doing three calculations, each with a boost in one of the Cartesian directions. However, alternatively one can also extract all data from a single calculation. As shown in the first remark, the directions of the transition dipoles can be evaluated along with the excitation energies from a single calculation without knowledge of the oscillator strengths. This is sufficient to calculate the angles ϕ0j be-

(1)

Note that the given accuracy of the energies and strengths does not represent the absolute accuracy of the TDDFT calculation with TDLDA but rather the accuracy of the evaluation scheme. A high accuracy is, e.g., important for the evaluation of transition densities as presented in the following.


10

Page 11 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For a specific transition 0 → n, eq (26) can be written as

tween boost and transition dipoles that appear in the inner products (k · r0j ) of eq (19). When this information is known, the missing oscillator strengths can be extracted from the line heights (ϑ) (f0j )γ . While this relieves one in principle of performing three calculations, it requires the boost direction ϑ to be non-orthogonal to any of the transition dipoles, and is probably less accurate. Third, we note that the scheme as presented here is designed to obtain dipole active excitations. However, the real-time approach can be extended to also identify dipole-forbidden states with zero oscillator strength. This can be achieved by using an excitation and recording an observable that break the dipole symmetry. An early example of how non-dipole active excitations can be obtained in real-time TDDFT was given in Ref. 28. Another example is given below in the Davydov-like evaluation of the coupling strength between a donor and an acceptor molecule, where we apply the boost in the donor half-space only and record the dipole moments of both half-spaces separately, cf. Section 5.

4

Im {δ˜ n (r, ω0n )} = −T (k · r0n ) ρ0n (r) (27) ∞ X sin ((ω0n − ω0j )T ) − (k · r0j ) ρ0j (r) . (ω0n − ω0j ) j=1,j6=n The appearance of the propagation time T in the first term of the right hand side of eq (27) ) → T for ω → 0. is due to sin(ωT ω The second term on the right hand side of eq (27) reveals how the transitions 0 → j with j 6= n show up in the Fourier transform that aims at showing the transition density for 0 → n. When |(ω0n − ω0j )T | ≫ 1 for all j 6= n, this term can be neglected so that ρ0n (r) ≈ − [T (k · r0n )]−1 Im {δ˜ n (r, ω0n )} , (28) which leads to the proportionality mentioned earlier in eq (12).

4.2

In the case of Na4 the transition densities of the two largest excitations Z1 and Y2 are shown in Figure 3 after a 50 fs propagation with a boost in the direction of the respective transition dipole. Both transition densities dominate Im {δ˜ n (r, ω)} at their respective excitation energy since they are strongly excited by the dipole-like boost (compare to Figure 1 (b) and (c)), i.e., we see that eq (28) holds. The left hand side of Figure 4, on the other hand, shows the transition densities approximated according to eq (28) for different propagation times and boost vectors for the Z2 excitation. We chose boost vectors in the yzplane such that excitations with a transition dipole perpendicular to the Na4 plane (i.e., X1, X2, and X3) are not excited. In the following, the term z-boost means k k eˆz while yz-boost ey + eˆz ). means k k √12 (ˆ From the spectrum in Figure 1 (b) from the 50 fs propagation it can be seen that especially the larger Z1 and, in the case of the yzboost, the Y2 excitation overlay the smaller Z2.

Evaluation of Transition Densities From Short Propagation Times

4.1

Accurate Transition Densities from Real-Time Propagation

Analytical Form of the Transition Densities obtained from a Real-Time Propagation

Not only the dipole moment, but also the timedependent density can be analyzed. If one replaces the dipole operator in eq (16) by the density operator n ˆ , one obtains the time-dependent density expressed as a sum over the transition densities. With a calculation that is similar to the one above and a spectral analysis under the same assumptions one obtains the Fourier transform of the density fluctuation (for ω > 0) Im {δ˜ n (r, ω)} = (26) ∞ X sin ((ω − ω0j )T ) − (k · r0j ) ρ0j (r) . (ω − ω0j ) j=1


11


Z1

Y2

Z2

Page 12 of 30

Boost ~z

a Figure 3: (a) Z1 and (b) Y2 transition densities from a 50 fs propagation after a boost in the direction of the respective transition dipole. (iso-value: 0.0001 a0 −3 ) Hence, the calculation of the transition density of the Z2 excitation is cumbersome, and the results obtained with different propagation times or boost directions differ very much, as seen in the left hand side column of panels in Figure 4. Due to the longest propagation time of 50 fs and z-boost, Figure 4 (a) is the most reliable of those on the left hand side. The transition density in Figure 4 (c) from a 10 fs propagation is strongly contaminated by the Z1 excitation, while the Y2 excitation is still switched off. Finally, in Figure 4 (e) also the Y2 excitation contributes and distorts the picture even more. The above issues can be addressed by a simple correction scheme that uses the accurate energies and strengths of the excitations from the improved evaluation scheme of Section 3.3. For a finite set of M close lying transitions 0 → n1...M , eq (27) can be written as a matrixvector equation with a M × M correlation matrix C and a diagonal normalization matrix A, i.e., (29) Im {δ n ˜} = C A ρ .

Boost ~yz

Figure 4: Z2 transition densities before (left column of panels) and after (right column of panels) the correction scheme from: (a) and (b) a 50 fs propagation with z-boost, (c) and (d) a 10 fs propagation with z-boost, and (e) and (f) a 10 fs propagation with yz-boost. (iso-value: 0.0001 a0 −3 ) previously, one simply has to invert the matrices C and A, i.e., ˜} . ρ = A−1 C −1 Im {δ n

(30)

The inversion of the diagonal matrix A is trivial and C can be inverted by standard algorithms. 62,82 Note that the approximation of eq (28) can be identified with setting C to the unit matrix. Note also that C depends only on the excitation energies while A depends only on their strengths and directions. This means that the knowledge of accurate excitation energies is already sufficient to obtain qualitatively correct

The vector Im {δ n ˜ } has as its entries the values of Im {δ˜ n} at energies ~ω0n1...M , while ρ holds the corresponding transition densities ρ0n1...M . A contains the normalization coeffi cients Ajj = −T k · r0nj . C contains the sin((ω0ni −ω0nj )T ) correlation factors Cij = (ω0n −ω and is 0nj )T i therefore symmetric with diagonal entries of 1. To evaluate the transition densities from the set of Im {δ˜ n} with much higher accuracy than


12

Page 13 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


results presented in this article, especially for the accurate calculation of coupling strengths presented in the next section. Appendix C addresses particularly difficult cases. The right hand side of Figure 4 contrasts the transition densities as obtained by using eq (30), including C, to the previously discussed ones obtained with the simple scheme. The transition densities from the new scheme are much more accurate and even with a short propagation time of 10 fs and a yz-boost that “mixes” all excitations, a reasonable transition density can be obtained, as seen in Figure 4 (f). To fully appreciate the strength of the new scheme one should realize that the Z2 excitation is hardly visible to the eye in the spectrum of Figure 1 (b) (which results from a 50 fs propagation). It is yet much less visible in a 10 fs propagation, but nevertheless the new scheme can extract the transition density reliably. For BChla the transition densities corresponding to the Q-band and CT transitions are shown in Figure 5 for a 50 fs propagation. For the calculation we chose the boost direction such that all of these states are excited, i.e., in the y direction of the coordinate system that the original structure 74 is given in. The Qy (a) and Qx (b) transitions are hardly affected by the correction. The uncorrected CT1 (c) and CT2 (e) transition densities, on the other hand, are highly contaminated, especially by Qy . After the correction both transition densities show less strength inside the bacteriochlorin ring, while CT1 (d) shows a higher strength at the top left and CT2 (f) a higher strength at the bottom right of the respective graphic. Summarizing this section we conclude that very often – as seen here for most of the excitations of Na4 and the spurious CT excitations of BChla – obtaining transition densities from real-time TDDFT reliably would require disproportionate, unreasonably long propagation times. The presented correction scheme solves this problem and returns, within the applied xc approximation, quantitatively accurate transition densities from the time-dependent density alone.

transition densities. With the following consideration we can devise a consistency check for the evaluation scheme that can also be turned into an additional step for yet further improving the accuracy of the transition densities. The transition dipoles used in the normalization matrix A stem from the evaluation of the dipole spectra discussed in Section 3.3 (in the following spec termed r0j ). On the other hand, one can also determine the transition dipoles after one has calculated the transition densities by using the equation Z dens r0j = rρ0j (r) d3 r . (31) Thus, the condition !

spec dens r0j = r0j

(32)

can be used as a consistency check: If the transpec dens sition dipoles r0j and r0j corresponding to a given excitation agree closely, then all steps spec of the scheme worked very well. If r0j and dens r0j have the same direction but differ slightly in their absolute values, i.e., Z spec dens r0j = Cr0j = C rρ0j (r) d3 r (33) with C ≈ 1, then the correction with C worked fine but the oscillator strengths (or |r0j |) used in the matrix A are not perfectly accurate. This can result from the finite accuracy with which spec r0j was obtained in the fit. Then, the transition densities and dipoles can be corrected for consistency by defining 1 spec consistent r0j = √ r0j C √ consistent ρ0j = Cρ0j ,

(34) (35)

which enforces the condition (32)(2)(3) . We always performed this consistency step for the (2)

corrected To see this, use r0j within A. Note that this can also be used to determine the oscillator strengths (or |r0j |) if they are not already known from the dipole moment analysis, which usually leads to C 6≈ 1. (3)


13


dipole coupling approximation, 46 which can also be translated into the real-time TDDFT framework. 49 The dipole-dipole coupling matrix element reads

50 fs

V dd = κ

(36)

(37)

κ is the orientation factor κ = eD · eA − 3(eD · eR )(eA · eR ). eA,D are the unit vectors in the direction of the respective transition dipoles of donor (D) and acceptor (A), and eR is the unit vector in the direction of R, which is the connection vector between the transition dipoles of donor and acceptor. ρA,D are the transition densities of the donor and acceptor excited states for which the coupling is to be evaluated. The second option is to go beyond the dipole coupling and use the full Coulomb coupling matrix element, which reads 46,47 ZZ e2 ρA (r) ρD (r) 3 3 ′ Coul V = d rd r . (38) 4πǫ0 |r − r′ |

d

f Figure 5: Transition densities of the Q-band and spurious CT transition of BChla from a 50 fs propagation. The Qy (a) and Qx (b) transition densities are hardly affected by the correction. For the CT states, uncorrected (c) and (e) and corrected (d) and (f) ones are depicted. (iso-value: 0.0002 a0 −3 )

5.1

|µA | |µD | 4πǫ0 R3

with the transition dipoles Z µA,D = −erA,D = −e rρA,D (r) d3 r .

b

5

Page 14 of 30

The costly explicit evaluation of the double integral, which is also known as transition density cube (TDC) method, 83 can be sidestepped by calculating the energy of the first transition density in the Coulomb potential generated by the second one. The latter improves the accuracy with respect to the usual TDC method. Comparing eqs (36) and (38) one sees that the dipole-dipole approximation approximates a transition density by its transition dipole. The third option is using the Davydov splitting, 48 which in the real-time domain amounts to evaluating beat frequencies in the timedependent molecular dipole moments. 49 The Davydov splitting exploits that the excitations of two identical molecules couple pairwise and split up in the frequency domain. It therefore includes in a natural way coupling effects that go beyond the ones that are taken into account by the previous two options. In the present resonant case, i.e., donor and acceptor excitation

Real-Time Evaluation of Coupling Strengths Three Ways of Evaluating Coupling Strengths

As a final application of the previously discussed concepts we calculate the energy transfer between molecules in real-time TDDFT. This is of practical relevance because energy transfer efficiencies between chromophores are decisive for understanding light-harvesting systems. In the following we evaluate coupling strengths in three different ways. First, one can use Förster’s famous dipole-


14

Page 15 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


lines get closer in frequency as the coupling strength between the chromophores becomes smaller (e.g. due to increasing separation). For very small couplings, the lines can hardly be distinguished from each other. However, with the high accuracy and spectral resolution that we can achieve with the above presented evaluation techniques, even small couplings can be evaluated from moderate propagation times, as we demonstrate with the following examples.

energies are the same (E ∗ ), the resulting excitations of the combined system have energies E ∗ ±V Dav so that the splitting equals two times the coupling strength V Dav . 49 If the two chromophores are aligned parallel to each other (as done later in this section), the two different energies that result from the splitting correspond to a symmetric and an antisymmetric coupling between donor and acceptor excitations. The antisymmetric one has a vanishing transition dipole such that the corresponding line cannot be seen in the spectrum of the combined system. This is derived and further discussed by means of an exciton model Hamiltonian in Appendix D. Since the energies of the combined system lie symmetrically with respect to the excitation energy of a single uncoupled system (E ∗ ± V Dav ), one can in principle evaluate the coupling strength between two identical molecules in the following way: First calculate the excitations of a single molecule, then calculate the excitations for the two coupled molecules (of which only the symmetrically coupled ones are visible in the dipole signal), and finally compare the energy of the symmetric excitation of the coupled system with the excitation energy of the single molecule. The difference equals V Dav . However, in practice we found this method to be prone to inaccuracies as small differences in the numerical representation of the single-molecule and the two-molecule system can heavily affect the calculated coupling strengths (also see Appendix E). Instead it is favorable to extract the coupling strength from a single calculation of the combined system and compare the energies of the symmetric and antisymmetric excitations. Still, as already mentioned, the antisymmetric line can neither be seen in the total spectrum nor be excited by a global dipole-like excitation. This problem can be circumvented, however, by applying the boost only to one of the molecules(4) and by integrating the dipole moment only over the acceptor or donor half-space. Then both lines are visible, cf. Appendix D. The two

5.2

Test 1: Sodium Dimers

We first investigate two Na2 molecules that are both aligned with their symmetry axis in z direction and separated in x direction such that their single transition dipoles are perpendicular to the connection vector R = Rex . Na2 is a good model for a two-level system since it shows one strong excitation along the symmetry axis of the dimer, 76 and therefore is a natural, transparent first test case. Figure 6 (a) shows the coupling strengths between the sodium dimers using the above described approximations for different distances R on a logarithmic scale. The dipole-dipole data points show the expected R−3 behavior from eq (36). The transition dipoles used in the dipole-dipole coupling approach were calculated from the transition densities of the Na2 excitation according to eq (37). Hence, the coupling strengths from dipole-dipole and transition density coupling agree per construction for large distances(5) . Due to Na2 being close to a two-level system and the high accuracy of our real-time scheme, we could also evaluate the Davydov coupling strength from a T = 50 fs propagation even for distances up to 50 a0 , i.e., very small couplings (these date are not shown in Figure 6 for reasons of transparency). The energetic difference between the spectral lines in this case is only ≈ 0.003 eV while the line width is Tπ ≈ 0.04 eV and thus larger by an order of magnitude. For distances of less than 20 a0 the dipoledipole coupling begins to differ significantly from the full Coulomb coupling. Figure 7 shows

(4) In our case by applying the boost only in the negative x half-space (Na2 -System) or z half-space (BChla system) of our real-space grid.

(5)

We tested this explicitly for up to R = 1000 a0 to check for errors in the evaluation of eq (38).


15


Coul

)/V

dd

[ ] 0.05

10%

a

5%

10

11

12

14 16 Distance R [a0]

0.1

18

20

14

22

Dipole-Dipole (Qy-Qy) Full Coulomb (Qy-Qy) Davydov (Qy-Qy)

0.05

BChl-BChl

0.05

16 18 20

30 Distance R [a0]

40

50

60

Figure 7: Error of the dipole-dipole approximation compared to the full Coulomb coupling for the Na2 -Na2 and the BChla-BChla donoracceptor pairs.

b 0.01

50%

Coul

Na2-Na2

[ ]

0.1

Na2 - Na2 BChl - BChl

100%

Dipole-Dipole Full Coulomb Davydov

(V -V

Coupling strength [eV]

0.5

Coupling strength [eV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 30

0.04

lar spectrum. This effect starts to set in at R = 11 a0 and is evident when one looks at the dipole spectrum (not shown here), because the spectral lines from the symmetric and the antisymmetric excitations start to have heavily unequal heights (cf. Appendix D) and further excitations appear nearby. The corresponding data points are therefore enclosed in brackets in Figure 6 (a). Given the three different evaluation schemes we can draw conclusions about the intermolecular separation at which the different coupling approximations are reliable: For the Na2 dimers, the dipole-dipole approximation reaches an accuracy of at least 10% for R > 20 a0 , pure Coulomb coupling of transition densities is reliable for R > 16 a0 , and the two-level approximation of the Davydovsplitting evaluation holds for R > 12 a0 .

0.03 0.001 0.02 12 10

12

14

16

14 16 18 20 30 Distance R [a0]

40

50

60

Figure 6: Coupling strengths between (a) two sodium dimers and (b) the Qy transitions of two BChla from dipole-dipole, full Coulomb interaction and Davydov splitting as a function of the intermolecular distance. See main text for the explanation of data points in brackets. the relative error of the dipole-dipole coupling compared to the full Coulomb one. The error in this case follows a power law ∝ R−2 and amounts to 10% at R = 20 a0 . The coupling strengths from the Davydov splitting in Figure 6 (a) fit perfectly to those of the pure Coulomb coupling for R > 16 a0 (up to R = 50 a0 ). Since both evaluation schemes are quite different, this agreement underlines the accuracy of our numerical implementation and overall consistency. For R = 16 a0 the Coulomb coupling strength differs from the Davydov one by ≈ 3% and by ≈ 10% for R = 14 a0 . Finally, for R < 10 a0 , the evaluation based on the Davydov splitting was no longer possible as one cannot speak of separated dimers anymore, but rather gets a Na4 -type molecu-

5.3

Test 2: Bacteriochlorophylls

The question of energy transfer becomes of practical relevance as we turn to our second example, a system of two BChla aligned faceto-face to each other with nearly parallel Qy dipoles (κ = 0.975842) (cf. Appendix F for details of the setup and numerical details in general). The assumption of dipole-dipole coupling between BChla molecules has been made in the past, e.g., for setting up exciton models. 71 With the techniques presented above we can deter-


16

Page 17 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


distance of only about 15 a0 . At a closer distance of 11.5 a0 (cf. inset in Figure 6 (b)) the strengths predicted by the two methods differ by about 0.027 eV, i.e., by about 4.5% of the Davydov strength at this distance. One can hence conclude that the full Coulomb coupling mechanism for coupled systems of BChla is valid in a wide range of distances between the molecules, starting at about R > 11.5 a0 for this setup.

mine whether the dipole-coupling approximation is justified at a given intermolecular distance. The coupling strengths for distances between 14 a0 and 60 a0 from dipole-dipole and full Coulomb coupling are shown in Figure 6 (b). The total coupling strength for the BChla pair is smaller than for the Na2 pair, but qualitatively the situation is similar: Dipole-dipole coupling and full Coulomb coupling differ noticeably for small distances and get closer for large distances. For a quantitative statement about the range of validity of the dipole-dipole approximation, we again show its relative error in Figure 7. One finds the same R−2 decay as for the Na2 -Na2 system, but the absolute error is much larger for the BChla-BChla system. Only for a separation as large as 50 a0 does the error fall below the 10% level which is indicated by the horizontal dashed line. This finding is plausible in view of the spatial extension of a BChla molecule. In addition one should note that the face-to-face setup that we chose is close to optimal for the dipole-dipole coupling. Therefore, R = 50 a0 can be seen as a lower bound for the distance at which the dipole-dipole approximation between BChla can be trusted. We expect that these conclusions would not change if the two molecules were embedded in an isotropic dielectric medium instead of vacuum, because the dielectric medium would scale the coupling but not change its character. Energy transfer models that rely on dipole-coupling between BChla at a closer separation than 50 a0 should therefore be viewed with caution. Finally, we also employed the Davydov splitting approach using the same setup at various distances R between 11 a0 and 31 a0 . The Davydov data points are marked with blue stars in Figure 6 (b). Due to numerical intricacies explained in Appendix E, the distances used for this approach are integer multiples of the grid spacing for an accurate evaluation. Details of the numerical setup and estimates of the expected accuracy are given in Appendix F. The coupling strength from the Davydov splitting and the full Coulomb coupling agree well for the larger distances, but also for a

6

Conclusion

In this article we showed how to evaluate electronic excitations and transition densities from real-time TDDFT with an accuracy that goes beyond the one reached with the typically used evaluation of the PSD or DSF. By deriving the analytical form that real-time TDDFT spectra take after a boost excitation and fitting these to the numerical simulation data, the accuracy can be considerably increased without increasing the simulation time and computational effort. We demonstrated that this allows to reveal within real-time TDDFT excitations that carry little oscillator strength and/or are energetically close to much stronger excitations. Furthermore we discussed the real-time simulation of energy transfer between two molecules, looking at the coupling between a model system of two Na2 dimers, and the more complex test case of two coupled BChla molecules. By evaluating the coupling in three different ways – dipole-dipole coupling approximation, full Coulomb coupling between two transition densities, and analyzing the Davydov splitting – we could determine the intermolecular separation at which each of these coupling schemes becomes (in)valid. In this way we demonstrated, among other things, that the dipole-dipole approximation for the coupling between two optimally aligned BChla is in error by more than 10% for intermolecular separations of less than 50 a0 . This result serves as a warning against the use of the dipole-dipole coupling approximation that is at the heart of Förster theory in BChla systems unless intermolecular separations are clearly larger than 50 a0 . On the other


17


hand we showed that the Coulomb coupling between two BChla from their transition densities is valid even for rather short distances and in perfect agreement with the one obtained from the Davydov splitting approach.

rough and the constraint therefore rather theoretical. In a less formal, physically motivated attempt for a tighter bound, one may argue that one can estimate hl|(k · rˆ)n |ji by C(kd)n , where C is a factor of order unity, such that kd ≪ 1 is a more useful condition. In practice, one can easily find the range of valid boost strengths, which usually is relatively wide, by varying the boost and testing the response for linearity. In the linear regime, the resulting energies and oscillator strengths are independent of the boost strengths. In the valid range, we generally chose a quite large boost to minimize perturbing influences from numerical noise, cf. Appendix F.

Appendix A

Approximations made for a Small Boost

The boost is applied in eq (13) by applying the operator eik·ˆr to the many-body groundstate wave function |0i. In our evaluation of observables the boost operator only appears in brackets of the form hl|eik·ˆr |ji = δlj + ik · rlj + O[hl|(k · rˆ)n≥2 |ji]. In eq (17) we neglect the latter terms and in eq (18) we neglect terms like (k · rlj )(k · rmn ) and higher ones with k · rlj = hl|(k · rˆ)|ji. In theory, one can estimate the above matrix elements by Z !n N X n ∗ 3N |hl|(k · rˆ) |ji| = k · ri ψ l ψ j d r i=1 ! n Z X N |k · ri | |ψl∗ ψj | d3N r ≤ i=1

≤ (kdN )

n

Z

|ψl∗ ψj | d3N r ,

Page 18 of 30

B

Details of the Spectral Analysis

The oscillator strengths appear as Fourier coefficients in the spectral analysis of eq (21). We use the Fourier transform in the definition Z F [g] = g˜(ω) = g(t)eiωt dt (40) Z 1 g˜(ω)e−iωt dω (41) F −1 [˜ g ] = g(t) = 2π and the fast Fourier transform algorithm 62 in our numerical calculations. The finite propagation time 0 → T can be expressed by a window function

(39)

WT (t) = Θ(t) − Θ(t − T ) .

where N is the number of electrons, k = |k|, d is the radius of a sphere that encloses the region of appreciable probability density of the finite system(6) , and ψl and ψj are the spatial representations of |li and |ji, respectively. In the first step of the estimate one uses the triangle inequality subsequently for the integral and the inner sum and in the second step the estimate k · ri ≤ kd. In this respect all neglected terms are of order O[(kdN )n≥2 ] and kdN ≪ 1 is a possible condition for a boost to be considered “small”. However, the estimate (39) is very

(42)

The convolution theorem states F Tr δµ(ϑ) (Θ(t) − Θ(t − T )) = (43) γ (ϑ) F Tr δµγ ∗ F [Θ(t) − Θ(t − T )] The Fourier transform of the sin-function is F [sin(ω0j t)] (ω) = i [δ(ω − ω0j ) − δ(ω + ω0j )] , 2

(44)

where δ(.) denotes Dirac’s delta distribution.

(6)

Note that one can always ensure that the system is centered around the origin of the coordinate system by a coordinate transformation r → r − r0 which must not r −r0 ) change the physics since eik·(ˆ only adds an additional global phase factor to the wave function.


18

Page 19 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the boost, such that it does not address Z1, e.g., by choosing the boost to be in the y direction. The stronger Y2 excitation is sufficiently higher in energy so that a 50 fs propagation is good enough to distinguish the excitations. The correction scheme works fine and does not change much, therefore (a) and (b) look very similar. However, if the initial boost is applied in yz direction, the Z1 transition density dominates the uncorrected transition density at the Y1 energy after a 25 fs propagation by far (Figure 8 (c)). In this case the correction scheme improves the picture a lot (Figure 1 (d)) but does not remove the Z1 perturbation completely from the Y1 transition density. Therefore, a somewhat larger propagation would be needed in this case.

The Fourier transform of the window is Z T eiωt dt F [WT ] = 0 ω sin T ω 2 = ei 2 T (45) ω 2 ω i sin ω T h ω 2 T + i sin T . = cos ω 2 2 2 The convolution of the Fourier transform of the dipole moment and the window is complex due to the finite propagation time T . Nevertheless, its imaginary part, i.e., the real cos-part of the window contribution in eq (45), is sufficient to determine the oscillator strengths. This part is also related to the dipole strength function, 17 as seen from eq (10). After simplification of the resulting product cos ω2 T sin ω2 T = 1 sin(ωT ), the convolution of both contribu2 shaped tions results in a spectrum with sin(ωt) ω peaks at the excitation energies. The spectrum for ω < 0 is simply the one for ω > 0 mirrored at ω = 0 and with negative signs. For |ω0j T | ≫ 1, i.e., propagation times much longer than the period of a dipole oscillation with frequency ω0j , the ω < 0-part, which provides no additional information, can be ignored. The graphical meaning of the upper condition ) is that the tails of the sin(ωT -functions for ω < 0 ω do not penetrate into the ω > 0 region.

C

Y1

Boost ~y

Boost ~yz

Evaluation of Difficult Transitions

d

As presented in Section 4 the method for evaluating transition densities works fine even for excitations such as the one that we labeled Z2 in the Na4 spectrum, although it is highly contaminated by nearby, larger excitations such as Z1 and Y2. However, we here point out that there are yet more difficult cases. In the following we consider the excitation that is labeled Y1 in the Na4 spectrum of Figure 1. It is completely hidden by the Z1 transition (if both transitions are addressed by an excitation). As can be seen in Figures 8 (a) and (b) the transition density of Y1 can easily be extracted if one chooses the initial excitation, in our case

Figure 8: Transition densities of Na4 before and after the correction: (a) and (b) from a 50fs propagation and y-boost as reference and (c) and (d) from a 25 fs propagation with yz-boost. (iso-value: 0.0001 a0 −3 )

D

Coupling Strengths Davydov Splitting

from

In this appendix we detail why evaluating the dipole moment of subsystems can reveal excitations that are not visible in the spectrum of the total system. One can easily understand the


19

Journal of Chemical Theory and Computation H-aggregate (V>0) DSF (z-component) [arb. units]

4

0

16 a0 z

-4

6 5

(a2)

02

01

x

Acceptor Fit

(b1)

2

-2

The Hamiltonian is expressed in the basis (|ADi, |AD∗ i, |A∗ Di), where |Di (|Ai) and |D∗ i (|A∗ i) are the ground state and excited state of the donor (acceptor). The aforementioned basis vectors are their antisymmetrized products. E0 is the ground-state energy, E ∗ is the energy of the excited state of the single donor and acceptor systems, and V is the coupling strengths between the excited states |AD∗ i and |A∗ Di. Furthermore, rD and rA are the transition dipoles of donor and acceptor, respectively. ˆ are E0 and E1,2 = E ∗ ± The eigenvalues of H V with corresponding eigenvectors |0i = |ADi and |1, 2i = √12 (|AD∗ i ± |A∗ Di). Note that |1i denotes the symmetric state, |2i the antisymmetric one. The following derivation assumes that donor and acceptor are separated from each other such that their states don’t overlap and exchange terms in the matrix elements below can be neglected. After a boost excitation of the donor (i.e. the boost is only applied in the donor halfspace) and under the same assumptions as in Section 3, the system is in the state

J-aggregate (V 0, left half of Figure 9) with a center-tocenter distance of 16 a0 , and once oriented inline (J-aggregate, V < 0, right half of Figure 9) with a center-to-center distance of 20 a0 (7) . In the sum of donor and acceptor spectra, the spectral line corresponding to the antisymmetrically coupled state at ω02 cancels out (especially for rA = rD ), which leads to the characteristic blue (red) shift in the absorption of H-aggregates (J-aggregates). Note that for the separate donor and acceptor dipole moments in eq (53), the two frequency components (ω01 and ω02 ) have equal weights. However, the line heights shown in each panel of Figure 9 are not equally n high, o because the (ϑ) (ϑ) figure does not show Im δ µ ˜γ , but Sγ ∝ n o (ϑ) ωIm δ µ ˜γ , i.e., there is an additional factor of ω. Note also that the evaluation of the acceptor dipole moment for the coupling strength is much more accurate for two reasons: First, one can easily distinguish close-lying lines on the acceptor side due to their unequal signs, which allows for an accurate fit. Second, since the donor half-space is excited, less perturbing states are visible on the acceptor side, namely only those that participate in energy transfer.

E

Figure 10: Representation of two identical molecules on the numerical grid. Details are explained in the text. coupling strength V . However, even if two molecules are chemically identical, their numerical representation on the real-space grid may not. An example is given in Figure 10. In (a) the two molecules (D and A) on the left and right hand sides of the grid are separated by 3.5 times the grid spacing ∆x. They are chemically identical but their representations on the grid, i.e., the relative positions of their atoms with respect to the grid points, are different, which also makes the two molecules numerically different. Most easily this can be seen if one calculates the excitation energies of a single molecule twice and shifts the atom coordinates in the second calculation by, e.g., 0.5∆x in some direction. The resulting excitation energies will numerically differ a little, which, carried over to the case of Figure 10 (a), leads to ∆E = 21 (ED∗ − EA∗ ) 6= 0. On the other hand, the two molecules in (b) are separated by 3∆x, i.e., an integer multiple of the grid spacing. In this case, their numerical representations on the grid are equal, which carries the analytically expected ∆E = 0 over to the numerical result. The data presented for the Davydov splitting between two BChla in Figure 6 (b) are all calculated at distances R that equal an integer multiple of the grid spacing to ensure the latter condition. At other distances between the molecules the coupling strengths resulting from the same evaluation (cf. Appendix D) may

Influence of the Molecular Representation on the Numerical Grid on the Davydov Splitting

The calculation of the coupling strengths from the Davydov splitting as presented in Section 5 and further explained in Appendix D, includes a numerical intricacy, especially for the BChla data presented in Figure 6 (b). As discussed above, the coupling between two identical twolevel molecules results in two spectral lines with an energetic separation of two times the (7)

(z)

In this special case Figure 9 (a1) and (a2) show Sz for the H-aggregate for donor and acceptor, respectively, (x) while (b1) and (b2) show Sx for the J-aggregate.


21


be significantly overestimated. The reason is that in the numerically non-resonant √ case with ∆E 6= 0 the line splitting equals 2 V 2 + ∆E 2 , which is larger than 2|V | in the analytically expected resonant case. E.g., shifting the BChla of Sections 3 to 5 by 0.5∆x in each direction results in ∆E ≈ 0.008 eV, which is small compared to the absolute accuracy of the simulation but significant compared to the coupling strengths reported in Figure 6. Still, if ∆E is known, e.g., from the calculation of the single molecule with the two different grid representations, this situation can be remedied.

F

Page 22 of 30

of the results at a distance of 20 a0 for both grid sizes. For the propagation we applied a boost of 0.0001 Ry in the negative x half-space (donor) in z direction (along the Na2 axis). Further propagation parameters were chosen as for a single Na2 . For calculating the BChla spectrum and transition densities we used a grid spacing of 0.18 ˚ A(≈ 0.34 a0 ). For a single BChla we used a grid radius of 22 a0 . We propagated for 50 fs with a boost strength of 0.01 Ry. Reducing the boost strength to 0.000001 Ry results in a deviation of 0.8% of the oscillator strength of Qy and of 0.07% of its excitation energy but yields line shapes with less quality due to the higher ratio of numerical noise. In order to calculate the coupling strengths between two BChla from the Qy transition density, we separated this transition density and its duplicate along the vector (1, 1, −1), which results in a face-toface orientation, with the distances reported in Figure 6 (b). Variations of numerical parameters such as rotating or shifting the molecule with respect to the grid or increasing the grid size by 4 a0 had an overall effect of reducing the resulting coupling strengths by only 0 − 2%. For the calculation of Davydov splittings with two BChla we rotated the original BChla coordinates 74 by +45◦ around the z axis and subsequently by −45◦ around the x axis. This moves the bacteriochlorin ring of a single BChla into the xy plane and the Qy transition dipole approximately along (1, −1, 0). Separating the two BChla along the z direction results in exactly the same face-to-face orientation of the two chromophores as used when calculating the coupling strengths from the transition densities above. For consistency, we chose, if possible, the same numerical parameters as used for the transition density calculation above. For the combined system we used a grid size of 39 a0 along the intermolecular connection vector (z axis) and 32 a0 in perpendicular directions. Reducing the grid size by 4 a0 in each direction had an effect of only ≈ 0.5% on the resulting coupling strengths such that further increasing the grid size was assumed to have no further effect. The propagation time was 100 fs for each calculation and the boost was applied

Numerical Details

All calculations were done with a new pseudopotential real-space and real-time TDDFT code called BTDFT that we devised (cf. Appendix G). The norm-conserving TroullierMartins 84 pseudo potentials used throughout had cutoff radii of: C 1.09 a0 (s and p), H 1.39 a0 (s), O 1.10 a0 (s and p), N 0.99 a0 (s and p), Mg 2.18 a0 (s) and 2.56 a0 (p and d), and Na 3.09 a0 (s, p, and d). For the calculations on the Na4 cluster in Sections 3 and 4 we used a grid spacing of 0.9 a0 . The half-axes of the boundary ellipsoid were all 25 a0 . We used a time step of 0.01 fs and a boost strength of 0.0001 Ry. We found, also for other sodium clusters, that the grid size still has some effect on excitations above 3 eV for the TDLDA functional while the excitations at lower energies are very stable. The energies of those excitations above 3 eV still vary a little even for very large grids with a radius of up to 50 a0 . For the calculations for a single Na2 (oriented in z direction) for the transition densities in Section 5 we chose a grid spacing of 0.8 a0 , a grid size of 35 a0 along the x axis and 25 a0 in perpendicular directions, a propagation time of 50 fs, and other parameters as for Na4 above. For the Davydov calculations with two coupled Na2 separated in x direction we used a boundary radius of 25 a0 perpendicular to the intermolecular connection axis. Along the x direction we used a half axis of 35 a0 for distances of less than 20 a0 and a half-axis of 45 a0 for larger distances. We checked the consistency


22

with a strength of 0.01 Ry in the donor halfspace in (1, −1, 0) direction, i.e., approximately along the Qy transition dipole for a strong signal. The evaluation of the different components (x and y) of the resulting donor- and acceptor dipole moments for the Davydov splitting was consistent throughout and revealed variations of usually < 0.5%. The only exceptions are the data points at 16.20 ˚ A ≈ 30.6136 a0 , which have an accuracy of approximately ±4%, and at 12.06 ˚ A ≈ 22.7901 a0 with an error of ±2%. Due to the topic discussed in Appendix E, the coupling strengths shown in Figure 6 are calculated with BChla distances R that are a multiple of the grid spacing, specifically 6.12 ˚ A, ˚ ˚ ˚ ˚ ˚ 8.10 A, 10.08 A, 12.06 A, 14.04 A, and 16.20 A.

G

1 0.8

a

Padé approximation

1

100 fs 0.8

0.6 Oscillator strength [arbitrary units]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.4 0.2 0 1 0.8

Fit to analytical function

c

100 fs

0.6 Oscillator strength [arbitrary units]

Page 23 of 30

Padé approximation

b

50 fs CT1 missing

0.6 0.4

CT2 misshaped

0.2

0.4 0.2 0 1 0.8

Fit to analytical function

d

50 fs

0.6 0.4 0.2

0

0 1.5

1.6

1.7

1.8

1.9

2

Excitation energy [eV]

2.1

1.5

1.6

1.7

1.8

1.9

2

2.1


Figure 11: Dipole strength of BChla in Padé approximation (left half) and from our method (right half) with Lorentzian broadened lines from 100 fs (panels (a) and (c)) and 50 fs (panels (b) and (d)) propagations.

The BTDFT Program Set

We developed a new program set called BTDFT (Bayreuth DFT) for efficient, highly parallel real-time TDDFT calculations. BTDFT began as an extension to the real-space pseudopotential code PARSEC 85 and its real-time TDDFT extension. 28 It has by now been developed into an independent ground-state and real-time code with a much higher efficiency than its predecessors. Special emphasis has been put on optimizing the structure of the realspace grid and the memory access patterns for the solution of the real-time TDDFT equations. Furthermore, in addition to PARSEC’s grid parallelization an additional orbital parallelization was implemented into the propagator. To guarantee a consistent and stable propagation, each time step is done in a predictor-corrector used scheme, whereby the potential at t + ∆t 2 for the current corrector step is a linear interpolation between the potential at t and the one at t + ∆t from the previous predictor or corrector step, respectively. Self-consistency is usually achieved with one or two corrector steps. This now allows us to address larger systems for longer propagation times, especially if the number of electrons involved is large. Like in PARSEC, core electrons are handled by norm-conserving Troullier-Martins pseudo potentials. 84 The real-space grid in BTDFT has an ellipsoidal shape. When zero-boundary con-

ditions apply, the boundary ellipsoid has to be large enough for the density of the finite system of interest to vanish in good approximation at the grid’s surface. For real-time TDDFT calculations that ionize the system absorbing boundary conditions can be used. 58

H

Comparison with the Padé approximant evaluation

As mentioned in the main text of the article, Bruner et al. 14 have recently presented a different approach for deducing information about electronic excitations from real-time propagations of limited length. Their method relies on equating the frequency dependent dipole moment with a Padé approximant. In this way, they can extract excitation energies from a real-time propagation with a given length with much greater accuracy than with the elementary evaluation of excitations described in subsection 2.3. Triggered by questions from a referee we here compare our approach to the one by Bruner et al. The left half of Figure 11 shows the frequencydependent dipole moment of the BChla molecule evaluated with the Padé approximant method of Bruner et al. for a propagation of 100 fs length (upper left panel (a)) and 50 fs


23

Journal of Chemical Theory and Computation 6

Refined DSF Fit with CT1

a

5

6 Dipole strength S(-hω) [1/eV]

(lower left panel (b)). For the evaluation we again used the trace of the dipole moments from different boost directions as explained in the main text so that the peak heights correspond to the oscillator strengths. As explained in Bruner et al., multiplying the time-dependent dipole signal with a decaying exponential function e−ηt before the Padé procedure prevents artifacts and leads to Lorentzian line shapes. In panels (a) and (b) of Figure 11 we used a decay rate of η = 0.001 [Rydberg time units]−1 (≈ 0.02 fs−1 ). The right half of Figure 11 shows the same dipole moment data, but evaluated with our approach as described in subsection 3.3, again for 100 fs (upper right panel (c)) and 50 fs (lower right panel (d)). For clarity, we broadened the excitations from our evaluation by Lorentzian functions with the same width as in the Padé evaluation. First comparing the two approaches for the long propagation run with 100 fs, i.e., panels (a) and (c), we see that both methods reveal the same excitations: The main excitations are found at 1.677 and 1.994 eV and the weak CT excitations are visible at 1.739 and 1.863 eV. However, now comparing the two methods for the shorter propagation time, panels (b) and (d), we see that the Padé approximant does not resolve the CT1 transition at 1.739 eV, but shows just one broad peak. Furthermore, it does not allow to accurately locate the CT2 transition at 1.863 eV. On the other hand, all of these excitations can unambigously be extracted from the same dipole data with a propagation time of 50 fs with our analysis, as seen in panel (d). Thus, our approach provides a better resolution. We note that this result does not change when a different line width (i.e. exponential decay rate in the Padé evaluation) is chosen. We also stress that our procedure of fitting the correct analytical form to the numerical data clearly shows when the fit is not complete. As a demonstration we show in Figure 12 the numerical dipole signal from the 50 fs propagation once with the complete fit with all four excitations (left panel (a)), and once the fit one obtains with the CT1 excitation missing (right panel (b)). Clearly, the fit is only accurate with the additional excitation included.

Dipole strength S(-hω) [1/eV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 30

50 fs

4 3 2 1 0

-1

Refined DSF Fit without CT1

b

5

50 fs

4

Fit without CT1 is not complete

3 2 1 0

-1 1.5

1.6

1.7

1.8

1.9

2


2.1

1.5

1.6

1.7

1.8

1.9

2

2.1


Figure 12: Fit to the frequency-dependent dipole signal of BChla with four excitations (panel (a)) and only three excitations (panel (b)). We believe that combining both methods – the Padé approximant and our analytical fit – may be the ideal approach for future real-time analysis, because the Padé approach allows to quickly determine the energies of the main excitations with little effort, and these can then serve as a first input for the accurate final analysis based on our approach. Acknowledgement We acknowledge support by the DFG program GRK1640 and the study program “Biological Physics” of the Elite Network of Bavaria.

References (1) Dreizler, R. M.; Gross, E. K. U. Density Functional Theory: An Approach to the Quantum Many-Body Problem; Springer: Berlin, Heidelberg, 1990. (2) Gross, E. K. U.; Dobson, J. F.; Petersilka, M. In Density Functional Theory II: Relativistic and Time Dependent Extensions; Nalewajski, R. F., Ed.; Springer: Berlin, Heidelberg, 1996; pp 81–172. (3) Marques, M.; Gross, E. Time-Dependent Density Functional Theory. Annu. Rev. Phys. Chem. 2004, 55, 427–455. (4) Marques, M. A. L.; Maitra, N. T.; Nogueira, F. M. S.; Gross, E. K. U.; Rubio, A. Fundamentals of Time-Dependent Density Functional Theory; Lecture Notes


24

Page 25 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in Physics; Springer: Berlin, Heidelberg, 2012; Vol. 837.

A. L.; Nogueira, F.; Oliveira, M. J. T.; Stewart, J. J. P.; Rubio, A. Insights into colour-tuning of chlorophyll optical response in green plants. Phys. Chem. Chem. Phys. 2015, 17, 26599–26606.

(5) Yabana, K.; Bertsch, G. F. Timedependent local-density approximation in real time. Phys. Rev. B 1996, 54, 4484– 4487.

(13) Li, Y.; Ullrich, C. A. The Particle-Hole Map: A Computational Tool to Visualize Electronic Excitations. J. Chem. Theory Comput. 2015, 11, 5838–5852.

(6) Tong, X.-M.; Chu, S.-I. Density-functional theory with optimized effective potential and self-interaction correction for ground states and autoionizing resonances. Phys. Rev. A 1997, 55, 3406–3416.

(14) Bruner, A.; LaMaster, D.; Lopata, K. Accelerated Broadband Spectra Using Transition Dipole Decomposition and Pade Approximants. J. Chem. Theory Comput. 2016, 3741–3750.

(7) Calvayrac, F.; Reinhard, P.-G.; Suraud, E.; Ullrich, C. Nonlinear electron dynamics in metal clusters. Phys. Rep. 2000, 337, 493–578.

(15) Theilhaber, J. Ab initio simulations of sodium using time-dependent densityfunctional theory. Phys. Rev. B 1992, 46, 12990–13003.

(8) Marques, M. A. L.; Castro, A.; Bertsch, G. F.; Rubio, A. octopus: a first-principles tool for excited electronion dynamics. Comput. Phys. Commun. 2003, 151, 60–78.

(16) Provorse, M. R.; Isborn, C. M. Electron dynamics with real-time time-dependent density functional theory. International Journal of Quantum Chemistry 2016, 116, 739 – 749.

(9) Marques, M. A. L.; Gross, E. K. U. In A Primer in Density Functional Theory; Fiolhais, C., Nogueira, F., Marques, M., Eds.; Lecture Notes in Physics; Springer: Berlin, Heidelberg, 2003; Vol. 620; pp 144–184.

(17) Yabana, K.; Bertsch, G. F. Timedependent local-density approximation in real time: Application to conjugated molecules. Int. J. Quantum Chem. 1999, 75, 55–66.

(10) Baer, R.; Siam, N. Real-time study of the adiabatic energy loss in an atomic collision with a metal cluster. J. Chem. Phys. 2004, 121, 6341–6345.

(18) Calvayrac, F.; Reinhard, P.; Suraud, E. Spectral Signals from Electronic Dynamics in Sodium Clusters. Ann. Phys. (N. Y). 1997, 255, 125–162.

(11) Andrade, X.; Strubbe, D.; De Giovannini, U.; Larsen, A. H.; Oliveira, M. J. T.; Alberdi-Rodriguez, J.; Varas, A.; Theophilou, I.; Helbig, N.; Verstraete, M. J.; Stella, L.; Nogueira, F.; Aspuru-Guzik, A.; Castro, A.; Marques, M. A. L.; Rubio, A. Real-space grids and the Octopus code as tools for the development of new simulation approaches for electronic systems. Phys. Chem. Chem. Phys. 2015, 17, 31371–31396.

(19) Castro, A.; Marques, M. A. L.; Rubio, A. Propagators for the time-dependent Kohn-Sham equations. J. Chem. Phys. 2004, 121, 3425–3433. (20) Chu, X.; Chu, S.-I. Role of the electronic structure and multielectron responses in ionization mechanisms of diatomic molecules in intense short-pulse lasers: An all-electron ab initio study. Phys. Rev. A 2004, 70, 061402–1 – 061402–4.

(12) Jornet-Somoza, J.; Alberdi-Rodriguez, J.; Milne, B. F.; Andrade, X.; Marques, M.


25


(21) Lein, M.; Kümmel, S. Exact TimeDependent Exchange-Correlation Potentials for Strong-Field Electron Dynamics. Phys. Rev. Lett. 2005, 94, 143003–1 – 143003–4.

Page 26 of 30

in real time. Phys. Rev. B 2007, 76, 035413. (29) Telnov, D. A.; Heslar, J. T.; Chu, S.I. Strong-field ionization of Li and Be: a time-dependent density functional theory with self-interaction correction. Chem. Phys. 2011, 391, 88–91.

(22) Thiele, M.; Gross, E. K. U.; Kümmel, S. Adiabatic Approximation in Nonperturbative Time-Dependent DensityFunctional Theory. Phys. Rev. Lett. 2008, 100, 153004–1 – 153004–4.

(30) Goncharov, V. A.; Varga, K. Real-space, real-time calculation of dynamic hyperpolarizabilities. J. Chem. Phys. 2012, 137, 094111–1 – 094111–7.

(23) Kammerlander, D.; Castro, A.; Marques, M. A. L. Optimization of the ionization time of an atom with tailored laser pulses: a theoretical study. Eur. Phys. J. B 2017, 90, 91–97.

(31) Wardlow, A.; Dundas, D. High-orderharmonic generation in benzene with linearly and circularly polarized laser pulses. Phys. Rev. A 2016, 93, 023428–1 – 023428–10.

(24) Vincendon, M.; Dinh, P. M.; Romaniello, P.; Reinhard, P.-G.; Suraud, É. Photoelectron spectra from full time dependent self-interaction correction. Eur. Phys. J. D 2013, 67, 97–103.

(32) Saalmann, U.; Schmidt, R. Non-adiabatic quantum molecular dynamics: basic formalism and case study. Zeitschrift für Phys. D Atoms, Mol. Clust. 1996, 38, 153–163.

(25) Dauth, M.; Graus, M.; Schelter, I.; Wießner, M.; Schöll, A.; Reinert, F.; Kümmel, S. Perpendicular Emission, Dichroism, and Energy Dependence in Angle-Resolved Photoemission: The Importance of The Final State. Phys. Rev. Lett. 2016, 117, 183001.

(33) Takimoto, Y.; Vila, F. D.; Rehr, J. J. Realtime time-dependent density functional theory approach for frequency-dependent nonlinear optical response in photonic molecules. J. Chem. Phys. 2007, 127, 154114–1 – 154114–10.

(26) De Giovannini, U.; Hübener, H.; Rubio, A. A First-Principles Time-Dependent Density Functional Theory Framework for Spin and Time-Resolved AngularResolved Photoelectron Spectroscopy in Periodic Systems. J. Chem. Theory Comput. 2017, 13, 265 –273.

(34) Wang, F.; Yam, C. Y.; Chen, G. Time-dependent density-functional theory/localized density matrix method for dynamic hyperpolarizability. J. Chem. Phys. 2007, 126, 244102–1 – 244102–10. (35) Meng, S.; Kaxiras, E. Real-time, local basis-set implementation of timedependent density functional theory for excited state dynamics simulations. J. Chem. Phys. 2008, 129, 054110–1 – 054110–12.

(27) Castro, A.; Appel, H.; Oliveira, M.; Rozzi, C. A.; Andrade, X.; Lorenzen, F.; Marques, M. A. L.; Gross, E. K. U.; Rubio, A. octopus: a tool for the application of time-dependent density functional theory. Phys. status solidi 2006, 243, 2465– 2488.

(36) Lopata, K.; Govind, N. Modeling Fast Electron Dynamics with Real-Time TimeDependent Density Functional Theory: Application to Small Molecules and Chromophores. J. Chem. Theory Comput. 2011, 7, 1344 – 1355.

(28) Mundt, M.; Kümmel, S. Photoelectron spectra of anionic sodium clusters from time-dependent density-functional theory


26

Page 27 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(37) Liang, W.; Chapman, C. T.; Li, X. Efficient first-principles electronic dynamics. J. Chem. Phys. 2011, 134, 184102–1 – 184102–7.

Robins, K. A.; Kirtman, B. Assessment of conventional density functional schemes for computing the polarizabilities and hyperpolarizabilities of conjugated oligomers: An ab initio investigation of polyacetylene chains. J. Chem. Phys. 1998, 109, 10489 – 10498.

(38) Ding, F.; Kuiken, B. E. V.; Eichinger, B. E.; Li, X. An efficient method for calculating dynamical hyperpolarizabilities using real-time time-dependent density functional theory. J. Chem. Phys. 2013, 138, 064104–1 – 064104–9.

(45) Kümmel, S.; Kronik, L. Hyperpolarizabilities of molecular chains: A real-space approach. Comput. Mater. Sci. 2006, 35, 321 – 326.

(39) Repisky, M.; Konecny, L.; Kadek, M.; Komorovsky, S.; Malkin, O. L.; Malkin, V. G.; Ruud, K. Excitation Energies from Real-Time Propagation of the Four-Component Dirac-Kohn-Sham Equation. J. Chem. Theory Comput. 2015, 11, 980 – 991.

(46) Förster, T. Zwischenmolekulare Energiewanderung und Fluoreszenz. Ann. Phys. 1948, 437, 55–75. (47) Sagvolden, E.; Furche, F.; Köhn, A. Förster Energy Transfer and Davydov Splittings in Time-Dependent Density Functional Theory: Lessons from 2Pyridone Dimer. J. Chem. Theory Comput. 2009, 5, 873–880.

(40) Dewhurst, J.; Krieger, K.; Sharma, S.; Gross, E. An efficient algorithm for time propagation as applied to linearized augmented plane wave method. Computer Physics Communications 2016, 209, 92– 95.

(48) Davydov, A. S. The Theory of Molecular Excitations. Sov. Phys. Uspekhi 1964, 7, 145–178. (49) Hofmann, D.; Körzdörfer, T.; Kümmel, S. Energy transfer and Förster’s dipole coupling approximation investigated in a realtime Kohn-Sham scheme. Phys. Rev. A 2010, 82, 012509.

(41) Casida, M. E. In Recent Advances in Density Functional Methods; Chong, D. P., Ed.; World Scientific: Singapore, 1995; pp 155–192. (42) Petersilka, M.; Gossmann, U. J.; Gross, E. K. U. Excitation Energies from TimeDependent Density-Functional Theory. Phys. Rev. Lett. 1996, 76, 1212–1215.

(50) Zangwill, A.; Soven, P. Density-functional approach to local-field effects in finite systems: Photoabsorption in the rare gases. Phys. Rev. A 1980, 21, 1561–1572.

(43) Casida, M. E.; Jamorski, C.; Casida, K. C.; Salahub, D. R. Molecular excitation energies to high-lying bound states from time-dependent densityfunctional response theory: Characterization and correction of the time-dependent local density approximation ionization threshold. J. Chem. Phys. 1998, 108, 4439–4449.

(51) Runge, E.; Gross, E. K. U. DensityFunctional Theory for Time-Dependent Systems. Phys. Rev. Lett. 1984, 52, 997– 1000. (52) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138.

(44) Champagne, B.; Stan, E. A. P.; van Gisbergen, J. A.; Baerends, E.-J.; Snijders, J. G.; Soubra-Gahoui, C.;

(53) Falke, S. M.; Rozzi, C. A.; Brida, D.; Maiuri, M.; Amato, M.; Sommer, E.; De


27


Page 28 of 30

Sio, A.; Rubio, A.; Cerullo, G.; Molinari, E.; Lienau, C. Coherent ultrafast charge transfer in an organic photovoltaic blend. Science 2014, 344, 1001–1005.

(63) Hofmann, D.; Körzdörfer, T.; Kümmel, S. Kohn-Sham Self-Interaction Correction in Real Time. Phys. Rev. Lett. 2012, 108, 146401.

(54) Perdew, J. P.; Wang, Y. Accurate and simple analytic representation of the electrongas correlation energy. Phys. Rev. B 1992, 45, 13244–13249.

(64) Vasiliev, I.; Ögüt, S.; Chelikowsky, J. R. Ab Initio Excitation Spectra and Collective Electronic Response in Atoms and Clusters. Phys. Rev. Lett. 1999, 82, 1919– 1922.

(55) Kümmel, S. Charge-Transfer Excitations: A Challenge for Time-Dependent Density Functional Theory That Has Been Met. Adv. Energy Mater. 2017, 7, 1700440.

(65) Andrae, K.; Reinhard, P.-G.; Suraud, E. Crossed Beam Pump and Probe Dynamics in Metal Clusters. Phys. Rev. Lett. 2004, 92, 173402.

(56) Schelter, I.; Foerster, J. M.; Ullmann, M.; de Queiroz, T. B.; Kümmel, S. unpublished

(66) Legrand, C.; Suraud, E.; Reinhard, P.-G. Comparison of self-interaction-corrections for metal clusters. J. Phys. B At. Mol. Opt. Phys. 2002, 35, 1115–1128.

(57) De Giovannini, U.; Varsano, D.; Marques, M. a. L.; Appel, H.; Gross, E. K. U.; Rubio, a. Ab initio angle- and energyresolved photoelectron spectroscopy with time-dependent density-functional theory. Phys. Rev. A 2012, 85, 062515.

(67) Mundt, M.; Kümmel, S.; van Leeuwen, R.; Reinhard, P.-G. Violation of the zero-force theorem in the time-dependent KriegerLi-Iafrate approximation. Phys. Rev. A 2007, 75, 050501.

(58) Dauth, M.; Kümmel, S. Predicting photoemission intensities and angular distributions with real-time density-functional theory. Phys. Rev. A 2016, 93, 022502.

(68) Wijewardane, H. O.; Ullrich, C. A. RealTime Electron Dynamics with ExactExchange Time-Dependent DensityFunctional Theory. Phys. Rev. Lett. 2008, 100, 056404.

(59) Li, X.; Smith, S. M.; Markevitch, A. N.; Romanov, D. A.; Levisbd, R. J.; Schlegel, H. B. A time-dependent Hartree-Fock approach for studying the electronic optical response of molecules in intense fields. Phys. Chem. Chem. Phys. 2005, 7, 233 –239.

(69) Thiele, M.; Kümmel, S. Photoabsorption spectra from adiabatically exact timedependent density-functional theory in real time. Phys. Chem. Chem. Phys. 2009, 11, 4631–4639.

(60) Mundt, M. Real-Time Approach to TimeDependent Density-Functional Theory in the Linear and Nonlinear Regime. J. Theor. Comput. Chem. 2009, 8, 561–574.

(70) Hofmann, D.; Kümmel, S. Self-interaction correction in a real-time Kohn-Sham scheme: Access to difficult excitations in time-dependent density functional theory. J. Chem. Phys. 2012, 137, 064117.

(61) Goings, J. J.; Lestrange, P. J.; Li, X. Realtime time-dependent electronic structure theory. WIREs Comput. Mol. Sci. 2018, 8, 8:e1341.

(71) Hu, X.; Ritz, T.; Damjanovic, A.; Autenrieth, F.; Schulten, K. Photosynthetic apparatus of purple bacteria. Q. Rev. Biophys. 2002, 35, 1–62.

(62) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. In Numerical Recipes in Fortran, 2nd ed.; Press, W. H., Ed.; Cambridge University Press, 1992.


28

Page 29 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(72) Cogdell, R. J.; Gall, A.; Köhler, J. The architecture and function of the lightharvesting apparatus of purple bacteria: from single molecules to in vivo membranes. Q. Rev. Biophys. 2006, 39, 227– 324.

(81) Dreuw, A.; Head-Gordon, M. Failure of Time-Dependent Density Functional Theory for Long-Range Charge-Transfer Excited States: The ZincbacteriochlorinBacteriochlorin and BacteriochlorophyllSpheroidene Complexes. J. Am. Chem. Soc. 2004, 126, 4007–4016.

(73) Cheng, Y.-C.; Fleming, G. R. Dynamics of Light Harvesting in Photosynthesis. Annu. Rev. Phys. Chem. 2009, 60, 241– 262.

(82) Anderson, E.; Bai, Z.; Bischof, C.; Blackford, S.; Demmel, J.; Dongarra, J.; Du Croz, J.; Greenbaum, A.; Hammarling, S.; McKenney, A.; Sorensen, D. LAPACK Users’ Guide, 3rd ed.; Society for Industrial and Applied Mathematics: Philadelphia, PA, 1999.

(74) Oviedo, M. B.; Negre, C. F. A.; Sánchez, C. G. Dynamical simulation of the optical response of photosynthetic pigments. Phys. Chem. Chem. Phys. 2010, 12, 6706–6711.

(83) Curutchet, C.; Mennucci, B. Quantum Chemical Studies of Light Harvesting. Chem. Rev. 2017, 117, 294–343.

(75) Broglia, R. A.; Colò, G.; Onida, G.; Roman, H. E. In Solid State Physics of Finite Systems; Broglia, R. A., Ed.; Advanced Texts in Physics; Springer: Berlin, Heidelberg, 2004.

(84) Troullier, N.; Martins, J. L. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B 1991, 43, 1993–2006.

(76) Kümmel, S.; Andrae, K.; Reinhard, P.G. Collectivity in the optical response. of small metal clusters. Appl. Phys. B Lasers Opt. 2001, 73, 293–297.

(85) Kronik, L.; Makmal, A.; Tiago, M. L.; Alemany, M. M. G.; Jain, M.; Huang, X.; Saad, Y.; Chelikowsky, J. R. PARSEC - the pseudopotential algorithm for realspace electronic structure calculations: recent advances and novel applications to nano-structures. Phys. status solidi 2006, 243, 1063–1079.

(77) Mundt, M.; Kümmel, S. Optimized effective potential in real time: Problems and prospects in time-dependent densityfunctional theory. Phys. Rev. A 2006, 74, 022511. (78) Williams, T.; Kelley, C. Gnuplot 4.6: an interactive plotting program. http:// www.gnuplot.info/, 2013. (79) Frigaard, N.-U.; Larsen, K. L.; Cox, R. P. Spectrochromatography of photosynthetic pigments as a fingerprinting technique for microbial phototrophs. FEMS Microbiol. Ecol. 1996, 20, 69–77. (80) Tozer, D. J. Relationship between longrange charge-transfer excitation energy error and integer discontinuity in KohnSham theory. J. Chem. Phys. 2003, 119, 12697–12699.


29


for Table of Contents use only Title: Access to Challenging Electron Dynamics by Accurate Evaluation of Real-Time Density Functional Theory Calculations Authors: Ingo Schelter ([email protected]) Stephan Kümmel ([email protected]) ToC graphic:


30

Page 30 of 30

Accurate Evaluation of Real-Time Density Functional Theory Providing

Recommend Documents