Performance of Tamm-Dancoff approximation on nonadiabatic couplings by timedependent density functional theory Chunping Hu, Osamu Sugino, and Kazuyuki Watanabe Citation: The Journal of Chemical Physics 140, 054106 (2014); doi: 10.1063/1.4862904 View online: http://dx.doi.org/10.1063/1.4862904 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/5?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Valence excitation energies of alkenes, carbonyl compounds, and azabenzenes by time-dependent density functional theory: Linear response of the ground state compared to collinear and noncollinear spin-flip TDDFT with the Tamm-Dancoff approximation J. Chem. Phys. 138, 134111 (2013); 10.1063/1.4798402 Analytical Hessian of electronic excited states in time-dependent density functional theory with Tamm-Dancoff approximation J. Chem. Phys. 135, 014113 (2011); 10.1063/1.3605504 Nonadiabatic coupling vectors for excited states within time-dependent density functional theory in the Tamm–Dancoff approximation and beyond J. Chem. Phys. 133, 194104 (2010); 10.1063/1.3503765 Assessment of time-dependent density functional schemes for computing the oscillator strengths of benzene, phenol, aniline, and fluorobenzene J. Chem. Phys. 127, 084103 (2007); 10.1063/1.2761886 Excited state nuclear forces from the Tamm–Dancoff approximation to time-dependent density functional theory within the plane wave basis set framework J. Chem. Phys. 118, 3928 (2003); 10.1063/1.1540109

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.217.6.6 On: Sun, 31 Aug 2014 10:49:11

THE JOURNAL OF CHEMICAL PHYSICS 140, 054106 (2014)

Performance of Tamm-Dancoff approximation on nonadiabatic couplings by time-dependent density functional theory Chunping Hu,1,a) Osamu Sugino,2 and Kazuyuki Watanabe1 1 2

Department of Physics, Tokyo University of Science, 1-3 Kagurazaka, Shinjuku, Tokyo 162-8601, Japan Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan

(Received 9 August 2013; accepted 9 January 2014; published online 3 February 2014) The Tamm-Dancoff approximation (TDA), widely used in physics to decouple excitations and deexcitations, is well known to be good for the calculation of excitation energies but not for oscillator strengths. In particular, the sum rule is violated in the latter case. The same concern arises within the TDA in the calculation of nonadiabatic couplings (NACs) by time-dependent density functional theory (TDDFT), due to the similarities in the TDDFT formulations of NACs and oscillator strengths [C. Hu, H. Hirai, and O. Sugino, J. Chem. Phys. 127, 064103 (2007)]. In this study, we present a systematic evaluation of the performance of TDDFT/TDA for the calculation of NACs. In the cases we considered, including a variety of systems possessing Jahn-Teller and Renner-Teller intersections, as well as an example with accidental conical intersections, it is found that the TDDFT/TDA performs better than the full TDDFT, contrary to the conjecture that the TDA might cause the NAC results to deteriorate and violate the sum rule. The surprisingly good performance of the TDA for NACs is probably because the TDA can partially compensate for the local-density-approximation error and give better excitation energies in the vicinity of intersections of potential energy surfaces. Our study also shows that it is important to use the TDA based on the rigorous full-TDDFT formulation of NACs, instead of using it based on an alternative approximate formulation. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4862904] I. INTRODUCTION 1, 2

Density functional theory (DFT) has become a standard tool for the calculation of ground-state electronic structures and related properties. As the extension to the time domain, time-dependent DFT (TDDFT)3–9 has also achieved great success in studying the properties of excited states, and the methodologies of TDDFT are still undergoing rapid development. In the past decade, there have been continuous efforts to develop TDDFT methods for nonadiabatic quantum simulations. For this purpose, there are two quantities to be provided by TDDFT: adiabatic potential energy surfaces (PESs) and nonadiabatic couplings (NACs). To construct adiabatic PESs, excitation energies are required, and these can now be routinely and efficiently computed from the TDDFT linear response theory.5, 10 In comparison, the TDDFT methods for calculating NAC have been developed relatively recently.11–20 In the formulation initiated by Chernyak and Mukamel,11 and which was further developed by Hu et al. in the Casida formalism,13 NACs are rigorously calculated in a way similar to that for oscillator strength. The formulation is achieved by using the nuclear derivative of the Hamiltonian as the perturbation operator, instead of using the dipole operator, as is done in the case of oscillator strength. Such a formulation can be easily implemented, and the efficiency of TDDFT has been demonstrated in nonadiabatic quantum simulations of practical systems.14, 21, 22 a) Electronic mail: [email protected]

0021-9606/2014/140(5)/054106/9/$30.00

Although the above TDDFT method for NAC is rigorous in its theoretical foundation, its performance is inevitably affected by the approximations used in the practical calculations. One commonly used approximation is for the exchangecorrelation interaction, with the simplest form given by the local density approximation (LDA). Although the LDA has been reasonably good for numerous studies of molecules and solids at nearly equilibrium geometries, it has been found that for NAC, the performance of the LDA deteriorates when PES intersections are approached.12, 13 This is similar to the observation that the LDA cannot be reliably used to construct PESs around conical intersections.23 There is still much work to be done in this area, and, in particular, there is a need for more sophisticated functionals in order to attain better performance near the intersections. Another frequently used approximation is for the electron-nucleus interaction. Although in numerous studies, modern nonlocal pseudopotentials with high transferability have been shown to reproduce the all-electron results at a significantly reduced cost, they do not perform well when calculating NACs.15 The NACs calculated for different atoms of a molecule or cluster are not well balanced, and the sum of the NAC vectors often shows large deviations from the sum rule, which states that the sum of the NAC vectors on all the atoms should be zero, as required by translational invariance. In contrast, NACs can be accurately calculated in the all-electron atomic-orbital-basis framework.16 The pseudopotential problem for NAC has been analyzed by considering the locality of the perturbation and a way to avoid the problem has been developed: the nuclear derivative of the Hamiltonian

140, 054106-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.217.6.6 On: Sun, 31 Aug 2014 10:49:11

054106-2

Hu, Sugino, and Watanabe

∂H/∂R or the h operator is replaced by the nuclear derivative operator ∂/∂R or the d operator, and the d-matrix elements can then be evaluated with density functional perturbation theory.20 Bearing in mind the above background, in the present study we aim to evaluate another commonly used approximation in the TDDFT calculations, the Tamm-Dancoff approximation (TDA).24, 25 The TDA results in a Hermitian eigenvalue equation that can be solved using standard iterative techniques. The simpler form of the TDA equations leads to algorithms that have the advantages of better convergence, lower memory requirements, and only two (instead of three) superoperator-vector products per iteration.26–29 In the TDA, the quality of excitation energies is the same as in the full linear response calculations. Furthermore, the TDA can be used to circumvent the triplet instabilities.30–34 However, a wellknown disadvantage of the TDA is that the sum rules are no longer fulfilled, which leads to poor results in the calculated oscillator strengths.35–39 This then raises a critical issue: since NACs and oscillator strengths are calculated in similar ways, how accurately does the TDDFT/TDA calculate NAC? It is noted that previously Tavernelli et al.14, 17 and the Furche group18 have used the TDDFT/TDA for the calculation of NAC, but there has not been a systematic evaluation of the performance of the TDA near PES intersections. In this study, we compare the NAC results of the TDDFT/TDA with those of the full TDDFT for various Jahn-Teller and Renner-Teller systems. These systems possess intersections at known geometries and the properties of NACs in the vicinity of intersections are also well known, therefore, they are appropriate for benchmark calculations. We also study an example of accidental conical intersections for which the locations have to be numerically determined. The calculated results are somewhat surprising. In the cases we considered, we found that when intersection points were approached, the TDDFT/TDA gave more accurate results than did the full TDDFT; this was contrary to the conjecture that the results by the TDDFT/TDA would be poor, and the sum rule would be broken for NAC, as it is with oscillator strengths. In the literature, there have been similar reports that the TDDFT/TDA outperforms the full TDDFT for the calculation of excitation energies near conical intersections.31 The present results suggest that the TDA can be a good method to choose for the TDDFT calculation of NAC, and its improvement of the NAC results is because it partially compensates for the LDA error in the excitation energies without affecting other relevant quantities. Our study also shows the importance of using the TDA based on the rigorous full-TDDFT formulation of NAC, instead of using it based on an alternative approximate formulation. The present paper is organized as follows. In Sec. II, we introduce the TDDFT/TDA formulation for NAC, in both the h-matrix and d-matrix formulations. In Sec. III, implementation and computational details are given. In Sec. IV, we present the results of extensive calculations on molecular systems possessing various types of intersection, and compare the results of the TDDFT/ TDA and the full TDDFT. In Sec. V, we present our conclusions.

J. Chem. Phys. 140, 054106 (2014)

II. FORMULATION A. Full-TDDFT calculation of NAC from the Casida equation

The matrix element for NAC between the ground state and the excited states is defined by the following expression: τ0I,μ = 0 |

∂ |I , ∂Rμ

(1)

where  0 ( I ) is the many-body electronic wavefunction of the ground (Ith excited) state, and Rμ is the nuclear coordinate, in which μ represents the x, y, and z components and the atom index. In Ref. 13 we presented a formula for NAC, using the Casida formalism of TDDFT5 as follows: τ0I,μ =

h†μ S−1/2 FI 3/2

ωI

(2)

,

where ωI is the excitation energy, Sij σ,klτ =

(fkτ

δσ,τ δi,k δj,l , − flτ )(εlτ − εkτ )

(3)

which is constructed from the occupation number fiσ and the Kohn-Sham (KS) eigenenergy εiσ , and hij σ,μ = ψiσ | hμ |ψj σ  = ψiσ |

∂H |ψj σ , ∂Rμ

(4)

which is a recast of the nuclear derivative of the many-body Hamiltonian in the KS density matrix. Here, the KS orbitals ψ iσ are assumed to be real. FI is the renormalized eigenvector of the Casida equation,5, 10 i.e.,

FI = ωI2 FI ,

(5)

where



ij σ,klτ = δσ,τ δi,k δj,l (εlτ − εkτ )2 + 2 (fiσ −fj σ )(εj σ −εiσ )  × Kij σ,klτ (fkτ − flτ )(εlτ − εkτ ), (6)

with K being the KS matrix of the Hartree and exchangecorrelation kernel. The renormalization condition requires † that FI FI = 1 if is independent of ωI . The formulation of NAC is very similar to that of oscillator strength, and requires only that the dipole operator is replaced by the h operator for the perturbation. Based on the above h-matrix formulation, further derivation using the d-matrix19 gives NAC as τ0I,μ = ωI d†μ S1/2 FI , 1/2

(7)

where dij σ,μ = ψiσ | dμ |ψj σ  = ψiσ |

∂ |ψj σ . ∂Rμ

(8)

The d-matrix formulation is derived directly from the original h-matrix formulation by using the relationship between the nuclear derivatives of the many-body and KS Hamiltonians. This avoids the problem of the pseudopotential approximation, and thus allows accurate results in the pseudopotential framework. The d-matrix can be calculated analytically using density functional perturbation theory.20

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.217.6.6 On: Sun, 31 Aug 2014 10:49:11

054106-3

Hu, Sugino, and Watanabe

J. Chem. Phys. 140, 054106 (2014)

B. Full TDDFT calculation of NAC in the iterative scheme

So far, our formula for NAC is based on the direct solution of the eigenvalue problem of the explicitly constructed

matrix [Eq. (5)]. On the other hand, the TDDFT calculations to obtain the solution to a low-lying excited state are usually done by using the following iterative method:40–42       A B X 1 0 X =ω , (9) B A Y 0 −1 Y where Aij σ,klτ = δi,k δj,l δσ,τ (εlτ − εkτ ) + Kij σ,klτ , Bij σ,klτ = Kij σ,lkτ .

(10)

(11)

The difference between the occupation number of occupied and virtual orbitals is assumed to be 1. In each iteration, it is only necessary to calculate the product of a matrix and a vector, not the matrix itself, and this is done by using the following identities: (A + B) |X + Y = ω |X − Y

(12)

(A − B) |X − Y = ω |X + Y .

(13)

and

Using the relationships S = (A − B)−1

(14)

FI = ωI (A − B)−1/2 |X + YI ,

(15)

and 1/2

the h-matrix formula of NAC in the iterative scheme of TDDFT can be expressed as τ0I,μ =

h†μ

|X + YI ωI

τ0I,μ =

,

(16)

X + Y|X − YI = 1.

(17)

This procedure for calculating NAC remains similar to that for calculating the oscillator strength. Equivalently, the d-matrix formulation of NAC in the iterative scheme gives τ0I,μ = d†μ |X − YI .

(18)

C. TDDFT formulation of NAC within the TDA

As a frequently used approximation to TDDFT, the TDA neglects the B matrix in Eq. (9), resulting Y = 0,43 thus the eigenvalue equation is simplified to (19)

Unlike the full-TDDFT equation, which gives excitation energies of ±ω, the above TDDFT/TDA equation gives only

h†μ |XI ωI

(20)

,

and the d-matrix formulation of NAC becomes τ0I,μ = d†μ |XI ,

(21)

where the renormalization condition is X|XI = 1. D. Approximate TDDFT formulation of NAC within the TDA

As mentioned in our previous work,20 it is easy to become confused in the TDDFT formulation of NAC, by the evaluation of the nuclear derivative of the many-body Hamiltonian in the KS system. One might be tempted to regard it as identical to the nuclear derivative of the KS Hamiltonian, since the many-body system is recast into the KS framework in DFT. We have pointed out the logic problem with this in Ref. 20. A recent work on NACs and Kohn-Sham derivative couplings, using a local orbital-basis set,44 has also paid attention to this aspect, where the many-body and Kohn-Sham Hamiltonians are explicitly treated as being different. On the other hand, the above logic was previously regarded as “possibly reasonable” if the TDA was used, but the practical performance of this has not been clarified. Therefore, we will analyze such an approximate TDDFT formulation of NAC within the TDA by the comparison with the “rigorous” TDDFT/TDA formulation. The replacement of hμ by hKS μ in Eq. (20) gives †

τ0I,μ =

where the renormalization condition is

AX = ωX.

positive values (+ω). This is because the physical image of the TDA is the decoupling of excitations and de-excitations, and this is usually good for the calculation of excitation energies and facilitates practical calculations. Within the TDDFT/TDA, the h-matrix formulation of NAC becomes

hKS μ |XI ωI

=

d†μ S−1 |XI ωI

.

(22)

This equation is referred to below as “h-by-hKS /TDA,” in order to distinguish it from the TDDFT/TDA equations in Sec. II C. Comparison of Eqs. (21) and (22) reveals that the difference comes from the powers of S and ωI . Another comparison of Eqs. (7) and (21) shows that in the TDDFT/TDA, 1/2 ωI S1/2 → 1. In this regard, reasonable results might be obtained using the approximate formulation, Eq. (22), within the TDA. Nevertheless, considerable differences are expected to arise in molecular systems when the powers of S and ωI are important for the excited-state configuration. This will be shown in detail with the example of HC8 S in Sec. IV A. III. IMPLEMENTATION AND COMPUTATIONAL DETAILS

The implementation of the present TDDFT/TDA formulation of NAC is carried out using the d-matrix within density functional perturbation theory, on the basis of the ABINIT code,45 which is a planewave pseudopotential approach. Ideally, the TDDFT/TDA equation should be solved by iteration, however, for benchmark calculations, the existing codes for

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.217.6.6 On: Sun, 31 Aug 2014 10:49:11

054106-4

Hu, Sugino, and Watanabe

J. Chem. Phys. 140, 054106 (2014)

the diagonalization of the explicitly constructed matrix in Eq. (5) is used for the A matrix in Eq. (19). This is done because the computational advantages of TDA are well known, and the present study is focused on the accuracy of the TDA for the TDDFT calculation of NAC. All calculations are performed within the adiabatic LDA with spin polarization, using the Teter Pade parametrization.46 The Troullier-Martins pseudopotentials47 with nonlinear core correction,48 generated by Khein and Allan, are used for various atomic species. Only the point (k = 0) is taken into consideration in the k-point sampling, which corresponds to the use of real wavefunctions. Convergence parameters, such as the supercell size, number of unoccupied orbitals, and kinetic energy cutoff, are examined to ensure reasonably accurate results.

(a)

(b)

IV. RESULTS AND DISCUSSION

In this section, we present the results of calculations of various molecular systems possessing intersection points where the ground state and the first excited state of these molecular systems are degenerate. For convenience we first discuss the Renner-Teller systems and then the Jahn-Teller systems. Finally, we discuss an example of systems where the conical intersections have to be determined numerically. A. Renner-Teller systems

Figure 1 shows the atomic geometries of the RennerTeller systems that are used in the calculation of NACs. Renner-Teller systems are known to possess intersection points when the atomic geometry is linear. When one atom is rotated on a circular contour around the intersection point, the angular NACs are quantized to 1 provided that the contour radius is sufficiently small. This quantization behavior is well known for triatomic Renner-Teller systems, and it has

H

y

y

S

x

q

N

q

H

O

q

z (b)

y x C1

C8 H

z

O

been demonstrated to be present for polyatomic Renner-Teller systems.49 Bearing this in mind, we have chosen NH2 , HCCS, and HC8 S as typical Renner-Teller systems with increasing numbers of atoms. In all cases, the leftmost atom (labeled as atom 1) is rotated on a contour with radius q around the collinear z axis and the angular NAC can be computed as 0 |

C C H

O

z (a)

S

x

FIG. 2. Calculation results for the NACs in NH2 by the TDDFT/TDA (squares) and the full TDDFT (circles), with the atomic geometry given in Fig. 1(a) and with the contour radius decreased from 2.0 bohrs to 0.1 bohr. The angular NACs are shown in (a) and the absolute values of the sum of the x components of the NACs on the three atoms are shown in (b).

(c) FIG. 1. Geometry of the Renner-Teller system as one atom is rotated on a circular contour with radius q around the intersection point (labeled as O) on the z axis. Illustrations are given for (a) NH2 , (b) HCCS, and (c) HC8 S systems. The nuclear configuration at the intersection point is a linear atomic chain.

∂ |1  = qτ1x (φ = 0), ∂φ

(23)

where φ is the contour angle. We choose φ = 0 for calculating the angular NACs when atom 1 is on the y axis, since the angular NACs in Renner-Teller systems do not have φ dependence. Figure 2(a) shows the q dependence of the angular NAC in NH2 , which is set according to the geometry in Fig. 1(a). The internuclear distance along the z axis is rH-N = 1.95 bohrs. As q is decreased from 2.0 bohrs, the results from both the TDDFT/TDA and the full TDDFT show an increasing tendency toward the ideal angular NAC of 1.0, consistent with what has been reported in the literature for using the stateaveraged CASSCF calculations.50 The values of the angular NACs from both methods are close to each other. However, after q = 1.0 bohr, the full-TDDFT result shows a sudden decrease from 0.97 to 0.88, while the TDDFT/TDA result continues to approach 1.0 as q is further decreased to 0.5 bohrs,

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.217.6.6 On: Sun, 31 Aug 2014 10:49:11

054106-5

Hu, Sugino, and Watanabe

J. Chem. Phys. 140, 054106 (2014)

and the result increases to 0.987. This value is nearly maintained at q = 0.2 bohrs. The performance of the full-TDDFT calculations is expected to be degraded when q becomes small, because we used the adiabatic LDA, which is known to give deteriorated results near intersection points.12, 13 Nevertheless, when the TDA is used, the performance is shown to be better than that of the full TDDFT, since more accurate NACs were obtained near to the intersection points. Figure 2(b) shows the sum (in absolute value) of the x component of the NACs for the three atoms of NH2 ; this can be used to judge the accuracy of the NAC results by considering the deviation from the ideal value of 0, as required by the sum rule. On the whole, the TDDFT/TDA and the full TDDFT give results at the same level, while the former better satisfies the sum rule as q ranges from 0.5 to 2.0 bohrs. Compared with the TDDFT calculation of oscillator strength, where the sum rule is violated when the TDA is applied, the present TDDFT/TDA calculation shows that such a problem does not arise for NAC, although the NACs were originally derived in a way similar to that for oscillator strength. Next we examine a four-atom Renner-Teller system, HCCS.51, 52 The atomic geometry is set according to Fig. 1(b) with rH-C = 2.0 bohrs, rC-C = 2.4 bohrs, and rC-S = 3.0 bohrs. Although there is no reference data for the direct computation of the angular NACs in HCCS, the quantization behavior of the angular NACs in the small q limit can be used for this, that is, the angular NAC should approach 1.0 as q decreases. Figure 3(a) shows the angular NACs in HCCS as a function of q. As q decreases from 2.0 bohrs to 0.5 bohrs, the results of both the TDDFT/TDA and the full TDDFT results increase, as can be seen in the case of NH2 . The TDDFT/TDA results

(a)

(b)

FIG. 3. Calculation results for NACs in HCCS by the TDDFT/TDA (squares) and the full TDDFT (circles), with the atomic geometry given in Fig. 1(b) when the contour radius q is decreased from 2.0 bohrs to 0.2 bohrs. The angular NACs are shown in (a), and the absolute values of the sum of the x components of the NACs on the three atoms are shown in (b).

TABLE I. The x components of the NACs, τ x (in bohr−1 ), for the individual atoms of HC8 S, calculated by the TDDFT/TDA, the full TDDFT, and the approximate h-by-hKS /TDA methods. The geometry of HC8 S is set as in Fig. 1(c), in which the carbon atoms are labeled in turn from left to right. The contour radius q is 1.0 bohr. Atom S C1 C2 C3 C4 C5 C6 C7 C8 H Sum

TDDFT/TDA

Full TDDFT

h-by-hKS /TDA

1.063 6.397 − 16.834 9.139 1.037 − 0.806 − 0.056 0.064 0.000 − 0.0025

1.678 10.097 − 26.572 14.426 1.636 − 1.273 − 0.088 0.102 0.000 − 0.004

0.303 2.737 − 6.956 3.865 0.435 − 0.332 − 0.023 0.029 0.002 − 0.0005

0.001

0.002

0.061

monotonically approach 1.0, up to q = 0.5 bohrs; beyond that, a positive deviation occurs. In contrast, the full-TDDFT results show that the deviations increase over the entire region. The deviation becomes several times larger than the quantized value of 1.0. This behavior is attributed to the deteriorating accuracy of the LDA. Nevertheless, to a large extent, the TDDFT/TDA improves the results of the full TDDFT and makes TDDFT more applicable in the smaller q region. Figure 3(b) shows the sum of the x components of the NACs for the four atoms of HCCS. The maximum value in the plot is less than 0.01 bohr−1 , which is regarded as a rather small value for practical applications. Therefore, the results of both the TDDFT/TDA and the full TDDFT satisfy the sum rule. To ensure the unexpectedly good performance of the TDDFT/TDA for the computation of the NAC, we further examined a ten-atom Renner-Teller system, HC8 S.52, 53 The atomic geometry is set according to Fig. 1(c), and the internuclear distance along the z axis is kept the same as it was in HCCS. The contour radius q is fixed at 1.0 bohr. The x components of the NACs, τ x , on individual atoms of HC8 S, are listed in Table I. Besides the results of the TDDFT/TDA and the full TDDFT, those from the approximate h-by-hKS /TDA method, using Eq. (22), are also listed for comparison. From τ x on atom S, the angular NACs are calculated to be 1.06 by the TDDFT/TDA, 1.68 by the full TDDFT, and 0.30 by the hby-hKS /TDA. Among these three methods, the TDDFT/TDA gives the result closest to the quantization value of 1.0. Also, the sum of τ x from the TDDFT/TDA is close to that of the full TDDFT, while the sum becomes much larger in the h-by-hKS / TDA. Further examination of the signs and the ordering of τ x for individual atoms in Table I shows that there might be a scaling relationship between the results from different methods. Figure 4 shows the ratio of τ x between the TDDFT/TDA and the full TDDFT, and that between the h-by-hKS /TDA and the full TDDFT. The atoms are labeled from 1 to 10 from the S atom to the H atom on the right. It can be clearly seen that the ratio of τ x is almost constant (∼0.63). Furthermore, the ratio of the excitation energy between the full TDDFT and the TDDFT/TDA is also ∼0.63, which is almost identical

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.217.6.6 On: Sun, 31 Aug 2014 10:49:11

054106-6

Hu, Sugino, and Watanabe

J. Chem. Phys. 140, 054106 (2014)

y y x

r

r (a)

(b)

q

h h

to that of τ x . The meaning of this can be understood from a comparison of Eqs. (16) and (20), which are the respective h-matrix formulations of NAC in the full TDDFT and the TDDFT/TDA. The excitation energy is the denominator in these formulas; therefore, the same ratios as seen above reveal that the numerators can be quite consistently obtained by the full TDDFT and the TDDFT/TDA, and the final results for NACs are rescaled by the excitation energies. Since the TDA might improve the performance of the LDA for the computation of excitation energies near the intersection points,31 we can see how this might lead to an improvement in the accuracy of the NACs calculated by the TDDFT/TDA. In contrast, the approximate h-by-hKS /TDA shows quite a different behavior, as can be seen in Fig. 4. The ratios are inconsistent for different atoms, and they are very different from the ratio of the excitation energies. This suggests that a good performance can only be obtained from the TDA that is based on the rigorous TDDFT formulations.

B. Jahn-Teller systems

Figure 5 shows the atomic geometry of the Jahn-Teller systems that were in the calculation of NACs. The Jahn-Teller systems have conical intersections with high symmetry: an equilateral triangle for a three-atom system in (a), a centered square for a five-atom system in (b), and a square in the top plane in (c). When one atom is rotated around the intersection point on the circular contour with radius q and angle φ, the angular NACs are calculated from the x and y components of the NAC on the rotated atom as ∂ |1  ∂φ   ∂ ∂ |1  cos φ + 0 | |1  sin φ . = −q 0 | ∂x2 ∂y2

0 |

(24)

x

O

q

FIG. 4. Comparison of the NAC results in Table I, as evaluated by the ratio of τ x between the TDDFT/TDA and the full TDDFT, and that between the h-byhKS /TDA and the full TDDFT. The ratio of the excitation energies between the full TDDFT and the TDDFT/TDA is shown by the broken line.

O

q

r

O

r

(c) Na9: side view & top plane FIG. 5. Geometry of the Jahn-Teller system as one atom is rotated on a circular contour with radius q around the intersection point (labeled as O). Illustrations are shown for (a) Na3 , (b) Na5 , and (c) Na9 systems. The nuclear configuration at the intersection point is an equilateral triangle in (a), a centered square in (b), and a square in the top plane in (c).

Unlike Renner-Teller systems, the angular NACs in a triatomic Jahn-Teller system are quantized to 0.5 at the small-q limit, and its integration over the contour angle from 0 to 2π gives the geometric phase (a value of π ), which is important for many geometric phase effects.54 We first take Na3 as an example. The internuclear distance in Fig. 5(a) is rNa–Na = 6.0 bohrs. Setting φ = 0, we examine the q dependence of the angular NAC as q is decreased from 5.0 bohrs to 0.5 bohrs. The results are shown in Fig. 6(a). Previous calculation results by TDDFT within the modified linear response (MLR) scheme are also shown for comparison. It is noted that the TDDFT-MLR can improve the performance of the LDA near the intersection points, and it produces an angular NAC that approaches 0.5 as q is decreased. On the whole, the TDDFT/TDA results are in general agreement with those of the full TDDFT, and they are slightly closer to those of the TDDFT-MLR as q is decreased to less than 1.0 bohr. Figure 6(b) shows the absolute value of the sum of the x components of the NACs for three atoms of Na3 . Although the sum seems to increase with the TDDFT/TDA, the maximum value is around 0.01 bohr−1 , thus the level of accuracy is still reasonable. Next, we examine the Na5 system, which is in the geometry shown in Fig. 5(b) with rNa–Na = 6.0 bohrs. Besides the q dependence, we also examine the φ dependence of the angular NAC, as shown in Fig. 7. Our previous work49 has verified that the angular NACs in polyatomic Jahn-Teller systems do not converge to 0.5 and that they show oscillatory behavior no matter how small q is. This oscillatory behavior can be seen in Figs. 7(a) and 7(b). Although qualitatively similar, the results of the TDDFT/TDA are quantitatively different from those of the full TDDFT. Further examination of the absolute values

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.217.6.6 On: Sun, 31 Aug 2014 10:49:11

054106-7

Hu, Sugino, and Watanabe

J. Chem. Phys. 140, 054106 (2014)

(a)

(b)

FIG. 6. Calculation results for the NACs in Na3 by the TDDFT/TDA (squares) and the full TDDFT (circles), with the atomic geometry given in Fig. 5(a) when the contour radius is decreased from 5.0 bohrs to 0.5 bohrs. The contour angle φ is set to 0◦ . The angular NACs are shown in (a) and the absolute values of the sum of the x components of the NACs for the three atoms are shown in (b).

of the sum of the x components of the NACs for all atoms, as shown by Figs. 7(c) and 7(d), indicates that the results of both the full TDDFT and the TDDFT/TDA results satisfy the sum rule, since the maximum values in both plot are very small. To study the geometric phase (α), we integrate the angular NAC in Figs. 7(a) and 7(b) over φ from 0 to 2π . The results are shown in Table II. For q at 1.0 and 0.5 bohrs, the TDDFT/TDA gives α as 3.149 and 3.166, respectively, which are very close to the ideal value of π . In contrast, the α given by the full TDDFT deteriorates from 3.592 to 4.199 as q is decreased. This shows that the TDA can improve the performance of the full TDDFT in a systematic way, not just in a random way, as the conical intersection is approached. The physical reason behind this phenomena might be that the TDA improves the excitation energies near conical intersections by partially compensating for the LDA error.31 The above results show that the geometric phase in polyatomic Jahn-Teller systems is computed more accurately by the TDDFT/TDA than by the full TDDFT. As another example of this, we did a calculation similar to the above, for Na9 in the geometry shown in Fig. 5(c). It is noted that Na9 is set in a three-dimensional geometry and the top plane is chosen for the contour. To remove multiple degeneracies, the symmetry was lowered from the body-centered cubic by setting the distance h between the centered Na atom and the top/bottom plane to be 0.5 bohrs larger than it was in the body-centered cubic. The same rNa-Na = 6.0 bohrs in Na5 was used, while the contour radius was set at q = 0.75 bohrs. As shown in Fig. 8(a), the angular NACs calculated by the TDDFT/TDA and the full TDDFT show qualitatively similar but quantitatively different behaviors. Integration over the whole contour

FIG. 7. Calculation results for the NACs in Na5 by the TDDFT/TDA (squares) and the full TDDFT (circles), as a function of the contour angle φ. The atomic geometry is given in Fig. 5(b) and the contour radius q is set to either 0.5 or 1.0 bohr. The angular NACs are shown in (a) and (b), and the absolute values of the sum of the x components of the NACs for the five atoms are shown in (c) and (d).

gives the geometric phase α as 3.01 for the TDDFT/TDA and 3.95 for the full TDDFT. It is noted that α can be smaller than π if the conical intersection under study is influenced by other conical intersections;54, 55 however, the fact that α is considerably larger than π is just due to errors in the calculation. This TABLE II. The geometric phase α in Na5 for different contour radii q, calculated by the integration of the angular NAC in Fig. 7 over the contour angle from 0 to 2π . q 1.0 0.5

TDDFT/TDA

Full TDDFT

3.149 3.166

3.592 4.199

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.217.6.6 On: Sun, 31 Aug 2014 10:49:11

054106-8

Hu, Sugino, and Watanabe

J. Chem. Phys. 140, 054106 (2014)

Na

H q

y

R H

O x r

FIG. 9. The geometry of the NaH2 system as one H atom is moved on the contour around the intersection point (O), which corresponds to an isosceles triangle with C2v symmetry.

FIG. 8. Calculation results of the NACs for Na9 by the TDDFT/TDA (squares) and the full TDDFT (circles), as a function of the contour angle φ. The atomic geometry is given in Fig. 5(c), and the contour radius q is set to 0.75 bohr. The angular NACs are shown in (a), and the absolute values of the sum of the x components of the NACs are shown in (b).

shows again the improvement of the TDDFT/TDA over the full TDDFT in three-dimensional Jahn-Teller systems. Meanwhile, the TDDFT/TDA accurately satisfies the sum rule, as shown by the small absolute values of the sum of the x components of the NACs in Fig. 8(b).

fer. The contour integration gives the geometric phase α as 3.19 and 3.54 in the TDDFT/TDA and the full TDDFT, respectively. Compared with the ideal value of π , the present results show that the TDDFT/TDA also performs better in the vicinity of accidental conical intersections. Figure 10(b) shows the absolute values of the sum of the x components of the NACs. Although the sum is larger at several points for the TDDFT/TDA, the magnitude is still on the order of 0.01–0.02 bohr−1 ; thus, the level of accuracy is still reasonable.

C. Systems with accidental conical intersections

Besides the symmetry-required Jahn-Teller conical intersections, there are many accidental intersections in chemical systems.56, 57 The positions of such intersections can only be determined by numerical calculations. The identification of the conical intersections in these systems is important work and requires much labor. For this reason, we take a well-known system, NaH2 , for our study of NAC.58, 59 The “accidental” conical intersections in this system have been interpreted by Baer and co-workers to be the elliptic Jahn-Teller type,58 which was shown to be of C2v symmetry, as shown in Fig. 9. In our calculation, we take a fixed rHH value of 2.18 bohrs, and RNa–H2 is determined as 3.6127 bohrs, where the degeneracy of the donor and acceptor orbitals (with equal fractional occupation of 0.5) is achieved. In the TDDFT calculation, the C2v symmetry is broken and the degeneracy is lifted by moving the H atom on a circular contour around the intersection point with a radius q = 0.5 bohrs (with integer occupations used in the linear-response TDDFT calculations). The calculated angular NAC as a function of the contour angle φ is shown in Fig. 10. The positions of the peaks of the two curves are in agreement, but they quantitatively dif-

FIG. 10. Calculation results for the NAC in NaH2 by the TDDFT/TDA (squares) and the full TDDFT (circles), as a function of the contour angle φ. The atomic geometry is given in Fig. 9 and the contour radius q is set to 0.5 bohrs. The angular NACs are shown in (a) and the absolute values of the sum of the x components of the NACs are shown in (b).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.217.6.6 On: Sun, 31 Aug 2014 10:49:11

054106-9

Hu, Sugino, and Watanabe

V. CONCLUSION

We presented a systematic evaluation of the performance of the TDA for the TDDFT calculation of NACs between ground and excited states. In the cases we considered, the TDA performed better than the full TDDFT, contrary to the conjecture that the TDA might cause the NAC results to deteriorate and violate the sum rule. This was conjectured because it is well known that the sum rule of the oscillator strength is lost with the TDA, and the original formulation of NAC is similar to that for oscillator strength. The calculation results in the vicinity of Jahn-Teller and Renner-Teller intersections, as well as in an example of accidental conical intersections, show that the TDA can improve the accuracy of the NAC, since the TDA can give better excitation energies as a result of partially compensating for the LDA error when the intersections are approached. Our work also shows that the good performance of the TDA can only be achieved on the basis of the rigorous TDDFT formulation of NAC. We believe that validation of the good performance of the TDA facilitates future nonadiabatic quantum simulations with TDDFT.

ACKNOWLEDGMENTS

The authors thank Jun Haruyama for fruitful discussions. This work was supported in part by the Project of Materials Design through Computics: Complex Correlation and Non-Equilibrium Dynamics, a Grant-in-Aid for Scientific Research on Innovative Areas, and the Next Generation Super Computing Project, Nanoscience Program, MEXT, Japan. K.W. acknowledges partial financial support from MEXT through a Grant-in-Aid (No. 22104007). Parallel testing of our program was performed on the supercomputers of the Institute for Solid State Physics, University of Tokyo. 1 P.

Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 3 E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984). 4 K. Burke, J. Werschnik, and E. K. U. Gross, J. Chem. Phys. 123, 062206 (2005). 5 M. E. Casida, in Recent Advances in Density Functional Methods, Part I, edited by D. P. Chong (World Scientific, Singapore, 1995), p. 155. 6 K. Yabana and G. F. Bertsch, Phys. Rev. B 54, 4484 (1996). 7 I. Vasiliev, S. Ö˘ güt, and J. R. Chelikowsky, Phys. Rev. Lett. 82, 1919 (1999). 8 M. A. L. Marques, X. López, D. Varsano, A. Castro, and A. Rubio, Phys. Rev. Lett. 90, 258101 (2003). 9 A. Castro, M. A. L. Marques, and A. Rubio, J. Chem. Phys. 121, 3425 (2004). 10 C. Jamorski, M. E. Casida, and D. R. Salahub, J. Chem. Phys. 104, 5134 (1996). 11 V. Chernyak and S. Mukamel, J. Chem. Phys. 112, 3572 (2000). 12 R. Baer, Chem. Phys. Lett. 364, 75 (2002). 13 C. Hu, H. Hirai, and O. Sugino, J. Chem. Phys. 127, 064103 (2007). 14 E. Tapavicza, I. Tavernelli, and U. Rothlisberger, Phys. Rev. Lett. 98, 023001 (2007). 15 C. Hu, H. Hirai, and O. Sugino, J. Chem. Phys. 128, 154111 (2008). 16 C. Hu, O. Sugino, and Y. Tateyama, J. Chem. Phys. 131, 114101 (2009). 17 I. Tavernelli, E. Tapavicza, and U. Rothlisberger, J. Chem. Phys. 130, 124107 (2009). 18 R. Send and F. Furche, J. Chem. Phys. 132, 044107 (2010). 19 C. Hu, O. Sugino, H. Hirai, and Y. Tateyama, Phys. Rev. A 82, 062508 (2010). 2 W.

J. Chem. Phys. 140, 054106 (2014) 20 C.

Hu, T. Tsukagoshi, O. Sugino, and K. Watanabe, Phys. Rev. B 87, 035421 (2013). 21 H. Hirai and O. Sugino, Phys. Chem. Chem. Phys. 11, 4570 (2009). 22 E. Tapavicza, A. M. Meyer, and F. Furche, Phys. Chem. Chem. Phys. 13, 20986 (2011). 23 B. G. Levine, C. Ko, J. Quenneville, and T. J. Martinez, Mol. Phys. 104, 1039 (2006). 24 S. Hirata and M. Head-Gordon, Chem. Phys. Lett. 314, 291 (1999). 25 M. Huix-Rotllant, B. Natarajan, A. Ipatov, C. Muhavini Wawire, T. Deutsch, and M. E. Casida, Phys. Chem. Chem. Phys. 12, 12811 (2010). 26 J. Hutter, J. Chem. Phys. 118, 3928 (2003). 27 The computational advantage of the TDA depends on practical implementations. In an integral-direct state-of-the-art implementation, the TDA usually has no considerable advantage over the full TDDFT. This is discussed, e.g., by Weiss et al. in Ref. 28. 28 H. Weiss, R. Ahlrichs, and M. Häser, J. Chem. Phys. 99, 1262 (1993). 29 S. Grimme, J. Chem. Phys. 138, 244104 (2013). 30 M. E. Casida, F. Gutierrez, J. Guan, F.-X. Gadea, D. Salahub, and J.-P. Daudey, J. Chem. Phys. 113, 7062 (2000). 31 F. Cordova, L. J. Doriol, A. Ipatov, M. E. Casida, C. Filippi, and A. Vela, J. Chem. Phys. 127, 164111 (2007). 32 M. E. Casida, A. Ipatov, and F. Cordova, in Time-Dependent Density Functional Theory, edited by M. Marques, C. Ulrich, F. Nogueira, A. Rubio, and E. Gross (Springer, Berlin, 2006). 33 M. J. G. Peach and D. J. Tozer, J. Phys. Chem. A 116, 9783 (2012). 34 M. J. G. Peach, M. J. Williamson, and D. J. Tozer, J. Chem. Theory Comput. 7, 3578 (2011). 35 J. Oddershede and P. Jorgensen, J. Chem. Phys. 66, 1541 (1977). 36 T. D. Bouman and A. E. Hansen, J. Chem. Phys. 66, 3460 (1977). 37 F. Furche, J. Chem. Phys. 114, 5982 (2001). 38 M. Casida and M. Huix-Rotllant, Annu. Rev. Phys. Chem. 63, 287 (2012). 39 M. Grüning, A. Marini, and X. Gonze, Nano Lett. 9, 2820 (2009). 40 R. E. Stratmann, G. E. Scuseria, and M. J. Frisch, J. Chem. Phys. 109, 8218 (1998). 41 R. Bauernschmitt and R. Ahlrichs, Chem. Phys. Lett. 256, 454 (1996). 42 A. Dreuw and M. Head-Gordon, J. Am. Chem. Soc. 126, 4007 (2004). 43 F. Neese, Coord. Chem. Rev. 253, 526 (2009). 44 E. Abad, J. P. Lewis, V. Zobaˇ c, P. Hapala, P. Jelinek, and J. Ortega, J. Chem. Phys. 138, 154106 (2013). 45 X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, Ph. Ghosez, J.-Y. Raty, and D. C. Allan, Comput. Mater. Sci. 25, 478 (2002). The ABINIT code is a common project of the Université Catholique de Louvain, Corning Incorporated, the Université de Liège, Mitsubishi Chemical Corp., and other contributors (URL http://www.abinit.org). 46 S. Goedecker, M. Teter, and J. Hutter, Phys. Rev. B 54, 1703 (1996). 47 N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991). 48 S. G. Louie, S. Froyen, and M. L. Cohen, Phys. Rev. B 26, 1738 (1982). 49 C. Hu, R. Komakura, Z. Li, and K. Watanabe, Int. J. Quantum Chem. 113, 263 (2013). 50 G. J. Halász, Á. Vibók, R. Baer, and M. Baer, J. Chem. Phys. 124, 081106 (2006). 51 M. Desouter-Lecomte, D. Dehareng, B. Leyh-Nihant, M. Th. Praet, A. J. Lorquet, and J. C. Lorquet, J. Phys. Chem. 89, 214 (1985). 52 A. Mitrushchenkov, R. Linguerri, P. Rosmus, and J. P. Maier, Mol. Phys. 107, 1549 (2009). 53 S. Wang, J. Li, X. Guo, L. Jiang, and J. Zhang, J. Mol. Struct.: THEOCHEM 900, 118 (2009). 54 M. Baer, Beyond Born-Oppenheimer: Electronic Nonadiabatic Coupling Terms and Conical Intersections (Wiley, Hoboken/New Jersey, 2006). 55 G. Halász, Á. Vibók, A. M. Mebel, and M. Baer, Chem. Phys. Lett. 358, 163 (2002). 56 D. R. Yarkony, Rev. Mod. Phys. 68, 985 (1996). 57 For systems of N > 3 atoms, the distinction between symmetry-determined and accidental degeneracy becomes somewhat nebulous (as pointed out by Professor C. A. Mead). Starting at a degenerate symmetric point, there are many ways to choose coordinates that break the symmetry but preserve the degeneracy. So, looking at different parts of the same conical intersection, one can view it as either symmetry-determined or accidental. 58 Á. Vibók, G. J. Halász, T. Vèrte´si, S. Suhai, M. Baer, and J. P. Toennies, J. Chem. Phys. 119, 6588 (2003). 59 D. R. Yarkony, J. Chem. Phys. 84, 3206 (1986).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.217.6.6 On: Sun, 31 Aug 2014 10:49:11

Performance of Tamm-Dancoff approximation on nonadiabatic couplings by time-dependent density functional theory.

The Tamm-Dancoff approximation (TDA), widely used in physics to decouple excitations and de-excitations, is well known to be good for the calculation ...
977KB Sizes 1 Downloads 0 Views