電子励起の第一原理計算： 時間依存密度汎関数法から のアプローチ

物 理 学 オープンセミナー, 6 月 22 日 電 子 励 起 の 第 一 原 理 計 算 : 時 間 依 存 密 度 汎 関 数 法 から のアプローチ 渡 辺 研 助 教 胡 春 平

Part I. Background

第 一 原 理 計 算 とは 第 一 原 理 (first principles) 計 算 というのは もっとも 基 本 的 な 原 理 に 基 づく 計 算 という 意 味 で, 電 子 間, 原 子 核 間,および 電 子 原 子 核 間 のクーロン 相 互 作 用 から 出 発 し, 量 子 力 学 の 基 本 法 則 に 立 脚 した 電 子 状 態 理 論 を 使 って 電 子 分 布 を 決 め, 物 質 の 諸 性 質 を 計 算 することを 指 す 非 経 験 的 電 子 状 態 計 算 とも 呼 ばれる 理 学 のキーワード, http://www.s.u tokyo.ac.jp/ken/key/03_3.html

波 動 関 数 理 論 : 多 体 波 動 関 数 密 度 汎 関 数 法 (DFT) Ψ( r, r2 HΨ( r, r2,..., rn ) = EΨ( r1, r2,..., r 1 n,..., r 1 n 密 度 汎 関 数 法 (density functional theory, DFT):あらゆる 物 理 量 の 期 待 値 は 基 底 状 態 の 電 子 密 度 の 汎 関 数 (すなわち 関 数 の 関 数 )である 1 2 + v KS eff 2 KS KS ( r) ϕ ( r) = ε i i ϕi ( r) KS 2 ρ( r) = ϕi ( r) i ) ) Kohn Sham (KS) 方 程 式 多 体 問 題 を 有 効 場 ( effective potential, Veff ( r)) の 中 の 一 電 子 の 軌 道 運 動 ( r) に 関 する 問 題 に 厳 密 に 置 き 換 える ただし 電 子 が 避 け ϕ i 合 い 運 動 することから 電 子 交 換 相 関 効 果 を 別 途 考 慮 する 必 要 複 雑 な 多 体 問 題 を 高 精 度 かつ 効 率 的 に 扱 う Walter Kohn, 密 度 汎 関 数 法 の 開 発, 1998 年 ノーベル 賞 ( 物 理 学 者 だが 化 学 賞 を 受 賞 した)

電 子 励 起 DFTによる 基 底 状 態 の 計 算 : 物 質 の 電 子 状 態 安 定 構 造 動 力 学 熱 力 学 の 性 質 などの 解 析 および 予 測 に 威 力 を 発 揮 さまざまの 分 野 で 大 きな 成 功 を 収 めた しかし 従 来 のDFTは 基 底 状 態 の 理 論 であり 電 子 励 起 状 態 を 厳 密 に 扱 うために 拡 張 が 必 要 基 底 状 態 :エネルギーEの 最 小 値 を 持 つ 固 有 状 態 E E[ ρ( r)] の 極 値 を 変 分 法 で 決 める: = 0 ρ 励 起 状 態 :エネルギーが 高 い 固 有 状 態 ; 基 底 状 態 と 直 交 している ただし DFTで 基 底 状 態 の 多 体 波 動 関 数 を 求 めない; 既 知 量 は 基 底 状 態 の 電 子 密 度 DFTを 時 間 領 域 に 拡 張 する: 時 間 依 存 密 度 汎 関 数 法

Kohn Sham (1965) 時 間 依 存 密 度 汎 関 数 法 (Time Dependent DFT, TDDFT) DFT TDDFT Runge Gross (1984) KS 2 ρ( r) = ϕi ( r) i ρ() r υeff () r = υext() r + dr + υxc() r r r 1 2 KS KS KS + υeff () r ϕi () r = εi ϕi () r 2 υ xc () r = δ E xc [ ρ] δρ() r E xc : exchange correlation energy KS 2 ρ( r, t) = ϕi ( r, t) i ρ( r, t) υeff (,) r t = υext(,) r t + dr + υxc(,) r t r r KS 1 2 KS i ϕi (,) r t = + υeff (,) r t ϕi (,) r t t 2 υ xc δ Axc[ ρ] δexc[ ρt] (,) r t = δρ(,) r t δρ(,) r t A xc : exchange correlation action

TDDFTで 励 起 状 態 を 求 める TDDFT 線 形 応 答 理 論 DFT: υ ext () r ρ () r 1 対 1の 対 応 : 外 場 vs 電 子 密 度 TDDFT: υ ext (,) r t ρ (,) r t TDDFT LR: δυ ext (,) r t δρ (,) r t (LR: Linear response) The perturbed state is a mixed state of ground and excited state components 直 交 条 件 : Ψ Ψ = gr ex 0 δρ( ω) 1 対 1の 対 応 : 摂 動 vs 電 子 密 度 の 変 化 Gross-Kohn: PRL 55, 2850 (1985) 電 子 密 度 がどうのように 応 答 するかを 計 算 することで 励 起 状 態 が 求 められる δρ(ω)のpoleの 位 置 は 励 起 エネルギー

Realization of TDDFT-LR 実 時 間 (Real-time) Yabana-Bertsch scheme Yabana and Bertsch, Phys. Rev. B 54, 4484 (1996). 周 波 数 領 域 (frequency-domain) The Casida formalism M. E. Casida, in Recent Advances in Density Functional Methods, Part I, edited by D. P. Chong (Singapore, World Scientific, 1995), p. 155; Jamorski, Casida, and Salahub, J. Chem. Phys. 104, 5134 (1996). Common feature: Preparation of reference state (initial state) Reference-state wavefunction by static Kohn-Sham equations: { 1 2 v [ ]} ext v v 2 + [ ρ] + H[ ρ] + xc ρ ϕi = εϕ i i

Real-time TDDFT-LR Procedures Perturbation Apply a weak perturbation to be in linear response regime In u direction δ v ext Time propagation Response Propagate the wavefunction and keep track of the dipole moment α() t = ϕi() t ru ϕi() t i The Fourier transform of the dipole moment is the response we seek. The excitation energy can be determined from the position of the dipole response. α( ω) excitation energy

実 時 間 発 展 TDDFTのメリット Straightforwardなformulation (favored by physicists) 摂 動 を 直 接 に 与 えてその 応 答 の 時 間 変 化 を 見 るので 同 じ 手 順 で 強 い 外 場 を 与 えて 非 線 形 応 答 も 見 られる 分 子 動 力 学 (molecular dynamics, MD)と 組 み 合 わせ ることができる TDDFT-MD 渡 辺 研 春 山 (D1) 石 見 (M2) しかし 実 時 間 発 展 はかなり 時 間 ステップが 必 要 なので 計 算 時 間 がかかる

Fourier-domain TDDFT-LR: The Casida formalism: 固 有 値 方 程 式 この 固 有 値 方 程 式 を 導 出 するためには 時 間 発 展 + 摂 動 + 線 形 応 答 の 式 を 使 ったが 最 後 の 方 程 式 はこれを 一 切 使 わなくて 基 底 状 態 の 情 報 のみを 用 いる 時 間 発 展 をする 必 要 がないので 中 小 規 模 の 分 子 系 の 計 算 に 効 率 が 良 い

The Casida Formalism Time-dependent Kohn-Sham equations: 1 2 σ + veff (,) r t (,) (,) 2 ϕi σ r t = i ϕi σ r t t σ Z I σ veff (,) r t = + δ vappl(,) r t + vscf(,) r t I ri σ ρ( r, t) σ vscf(,) r t = dr + vxc ρ, ρ (,) t r r r δv (,) r t = δv (,) r t + δv (,) r t σ σ eff appl SCF フーリエ 変 換 : δ v ( ω) = δv ( ω) + δv ( ω) σ σ eff appl SCF

The linear response of density matrix to the applied field is f f δ P ω δ ω jσ iσ eff ijσ( ) = vij σ( ) ω ( εi σ ε jσ) (Time dependent perturbation theory) On the other hand, δ v ( ω) = K ( ω) δp ( ω) SCF ijσ ijσ, klτ klτ klτ f kτ f lτ 0 ω ( εk τ εl τ) appl στ, i, k j, l K ijσ, klτ Pkl τ = vij σ klτ fl τ f kτ δ δ δ ( ω) δ ( ω) δ ( ω)

Separation of i j and j i elements f kτ δ δ δ ε ε 2 ( ω) 2 δστ, δik, δ jl, appl ω (Re δpkl τ)( ω) = δυij σ ( ω) ( fk τ fl τ)( εk τ εl τ) flτ > 0 kτ lτ στ, ik, jl, Kij σ, klτ klτ fk τ fl τ S I Ω S P υ 2 appl ( ω ) (Re δ )( ω) = δ ( ω) S ijσ, klτ = δ δ δ στ, ik, jl, ( f f )( ε ε ) kτ lτ lτ kτ Pole of the response: 2 ω I Ω = 0

Casida s Equation 2 I ΩF = ω F I Ω = δ δ δ ( ε ε ) ijσ, klτ σ, τ i, k j, l iσ jσ I 2 The eigenvalue of Ω is the squared excitation energy! + 2 ( f f )( ε ε ) K ( f f )( ε ε ) iσ jσ jσ iσ ijσ, klτ kτ lτ lτ kτ

Application of the Casida formalism A standard approach to calculate spectra of molecules; Implemented in almost all quantum chemistry codes. Many publications on various topics: atoms, molecules, clusters, etc. Ref.: Vasiliev and Martin, PRA 69, 052508 (2004) TDLDA Exp. Absorption spectra of the sodium clusters calculated by TDDFT LR combined with various exchange correlation potentials.

Part II. 電 子 励 起 状 態 に 対 する 第 一 原 理 計 算 手 法 の 開 発 (1) TDDFT 修 正 線 形 応 答 理 論 (modified linear response) (2) TDDFTによる 非 断 熱 結 合 係 数 (NAC)の 計 算 手 法 (ポスドク 研 究 @ 東 大 物 性 研 & 物 質 材 料 研 究 機 構 )

Part IIの(1) TDDFT 修 正 線 形 応 答 理 論 (modified linear response theory) Why do we need to go beyond ordinary linear response? Basic idea of modified linear response Implementation in real time TDDFT LR Implementation in frequency domain TDDFT LR

TDDFT-LR: Success and challenges Extensively applicable to electronic excitations in atoms, molecules, clusters, etc, however Valence excitations Most successful Typically bonding to antibonding excitations among valence orbitals Good accuracy within 0.4 ev of experiments Sufficient to qualitatively identify the spectra of most intense transitions LDA, GGA Rydberg and charge transfer excitations Systematic underestimation of the excitation energies Up to several ev

Rydberg transitions Rydberg transition: An electronic transition described as promotion of an electron from a bonding orbital to a Rydberg orbital. Rydberg orbital: For an atom, an orbital with principal quantum number greater than that of any occupied orbital of the ground state. Ry Orbital energy Enl =, 2 ( n μ ) nl Ry: Rydberg constant μ : quantum defect For a molecular entity, a molecular orbital which correlates a Rydberg atomic orbital in an atomic fragment produced by dissociation. Typically the extension of the Rydberg orbital is large compared to the size of the atom or molecular entity. nl

Molecular orbitals of N 2 Rydberg orbitals are spatially extended! Continuum Ionization Potential 3pπ u Rydberg orbitals 3pπ u 3sσ g 3sσ g π g π g Valence orbitals σ g π u σ g π u σ u σ u σ g σ g

Conventional TDDFT-LR within LDA or GGA functionals has poor performance on predicting Rydberg transition energies Rydberg transition in N 2 TDDFT/B3LYP Expt. 3 1 1 3 Σ Σ Π Σ σ 3sσ + g g g σ 3sσ + g g g σ 3 pπ u g u σ 3 pσ + u g u 10.62 12.00 11.18 12.20 11.55 12.90 11.53 12.98 Y. Tawada et al., J. Chem. Phys. 120, 8425 (2004) Underestimation: 1.3 ev (Mean Absolute Error)

Charge-transfer (CT) excitations Spatially separated: example of diatomic molecule Underestimation of excitation energies from TDDFT-LR. Example: Failure of TDDFT to predict the excitation energies of charge-transfer states of C 2 H 4 -C 2 F 4 R = 4 Å Ground state CT excited state A. Dreuw, J. L. Weisman & M. Head Gordon, J. Chem. Phys. 119, 2943 (2003). Ethylene tetrafluoroethylene (C 2 H 4 C 2 F 4 ) dimer TDDFT/SVWN: the lowest CT state (HOMO of C 2 F 4 LUMO of C 2 H 4 ) : 1.71 ev Theoretical analysis: at least 9.54 ev Underestimation by about 7.8 ev!

Systematic underestimation of Rydberg and CT excitation energies by TDDFT LR is usually blamed on Inadequacy in the long range behavior of exchangecorrelation (XC) functional Remedies: Asymptotic correction of XC functional J. Chem. Phys. 109, 10180 (1998). J. Chem. Phys. 121, 655 (2004). Long range corrected exchange J. Chem. Phys. 120, 8425 (2004). Accurate exchange / XC J. Chem. Phys. 119, 6475 (2003). Mol. Phys. 103, 711 (2005). 1 v () xc r ( r ) r However, a different point of view may be necessary how to extract excitation energy from linear response?

Intuitive view Rydberg 遷 移 や 電 荷 移 動 励 起 は 電 子 状 態 の 大 きな 変 化 を 伴 う 励 起 である しかし 通 常 のTDDFT 線 形 応 答 理 論 は 基 底 状 態 だけを 注 目 している: 基 底 状 態 に 摂 動 を 与 えて その 線 形 応 答 をみる 摂 動 を 与 えることにより 基 底 状 態 と 励 起 状 態 の 混 合 状 態 が 得 られる しかし TDDFT ALDAではその 混 合 状 態 の 電 子 密 度 が 基 底 状 態 と 励 起 状 態 のmixing ratioに 依 存 している そ の 結 果 resonance frequency( 励 起 エネルギー)もその 混 合 比 に 依 存 している Our idea: 修 正 線 形 応 答 理 論 基 底 状 態 と 励 起 状 態 の 平 均 的 な 電 子 状 態 に 対 して 線 形 応 答 理 論 を 適 用 することにより 電 子 状 態 の 変 化 による 誤 差 を 最 小 限 に 抑 える 中 間 励 起 状 態 ( 励 起 した 電 子 数 ( 占 有 数 )が0~1である 状 態 )に 摂 動 を 与 え 線 形 応 答 から 励 起 エネルギーが 得 られ その 励 起 エネルギーの 平 均 値 (0から1の 積 分 あるいは0.5のmidpointの 値 )を 取 る

Analytical view Approximate Excitation Energy from TDDFT-LR M. Petersilka, U. J. Gossman & E.K.U. Gross, PRL 76, 1212, 1996: LR KS KS Ω ε j εk + Kkj + ϕk() rϕj() r fxc(, r r ) ϕk( r ) ϕj( r ) drdr Ω j k ALDA, LR K = 1 ϕ () r ϕ () r ϕ ( ) ϕ ( ) dd r r r r r r kj k j k j In adiabatic local density approximation (ALDA), f xc LDA δvxc[ ρ] (, rr ) = δ( r r ) fxc () r = δ( r r ) δρ(,) r t KS, LDA KS, LDA LDA ε j εk + K kj + ρ j() r fxc () r ρk() r dr

Comparison with EXX In exact exchange only (EXX) system, { KS EXX } { KS EXX ε } j j εk k EXX, LR,, Ω +Δ +Δ + K kj -J kj J kj = ρk() r ρ j( r ) dd r r r r Δ= ϕ v v ϕ ˆHF EXX i i x x i (Gonze and Scheffler, PRL 82, 4416, 1999) Ω ALDA,LR EXX,LR and Ω are in quite different form! ALDA, LR EXX Ω and Ω are quite different for Rydberg and charge-transfer excitations!

Comparison with SCF U On the other hand, the SCF method (total energy difference between the self consistent field calculations of the ground state and excited state) gives the excitation energy as i Ω ΔSCF, LDA KS, LDA LDA KS, LDA LDA { j j } { k k } LDA ρ ( r) f ( r) ρ ( r) dr J j = ε + U ε + U xc k 1 LDA Jii ρi( ) fxc ( ) ρi( ) d +, acceptor orbital 2 r r r r = 1 LDA Jii + ρi( ) fxc ( ) ρi( ) d, donor orbital 2 r r r r Self interaction energy (vanishing in SI free theories, such as Hartree Fock theory) kj

Modified linear response : Integrating over occupation number for the j th level from 0 to 1 and using Janak theorem, we get Ω ALDA, LR averaged on occupation number 1 = Ω 0 ALDA, LR df j occ ΔSCF, LDA =Ω + K kj + { KS, LDA LDA} { KS, LDA LDA ε } j U j εk Uk = + { KS, EXX } { KS, EXX ε } j j εk k LDA ρ ( r) f ( r) ρ ( r) dr j xc k + + EXX Ω = +Δ +Δ + K kj -J kj Explicitly similar expression: The difference between shifted Kohn Sham eigenvalues minus the electron hole Hartree- Fock interaction good performance of Modified LR within ALDA is expected! de df tot j occ = ε Janak j K theorem kj -J kj

Implementation of Modified LR In principle same procedures with ordinary LR Take average of response from intermediate excited states Some extensions might be necessary. Real time TDDFT Frequency domain TDDFT

MLR for Real-time TDDFT For different initial states (intermediate states connecting the ground and excited state), repeat the following procedure: perturbation real time propagation dipole response Necessary extension: orbital analysis of dipole response for the assignment of transitions

TDDFT orbital analysis Decompose orbital wavefunctions into the contribution of all orbitals ϕ () t = a n () t ϕ (0), a n () t = ϕ (0) ϕ () t d n m m m Dipole moment Define = mn m u n m m n n n* α () t = dmm am () t am () t fn n m, m ϕ (0) r ϕ (0) A t a t a t n n n* m() = n() m () d a t a t f + cc n m n n* mn m () n () n.. (Constant ) n Am n Comparing Am ( ω) at a certain ω can analyze the peaks in αω ( ). ( ω) n m transition: Transition energy, Transition probability

Implementation of Real time TDDFT Code: TM TDDFT MD First Principles Simulation tool for Electron Ion Dynamics Suzuki Trotter split operator propagator Plane wave basis + supercell Troullier Martins Norm conserving pseudo potentials with s nonlocality XC potential: Adiabatic LDA (Perdew Zunger) Yabana s perturbation scheme Sugino & Miyamoto, PRB 59, 2579 (1999); ibid, B66, 89901 (E) (2002).

Charge transfer excitation of HCl 3, 4 5: π σ * charge transfer from Cl to H 5 th 3 rd, 4 th 5 th 3 rd, 4 th Ground state Excited state

Excitation energy vs. fraction of excited electron 9.0 ω t (ev) 8.5 8.0 7.5 7.0 Smoothly increase 1 + 1 HCl: X Σ A Π transition Orbital 3, 4 5 6.5 ω a 0.0 0.2 0.4 0.6 0.8 1.0 ω A f q i i i i (three point): 7.74 ev = A i i fi (five point): 7.75 ev CI: 7.84 ev

Summarized results and comparisons of vertical excitation energies State Transition LR M-LR CI SCF Exp. 1 a Π g V : σ g π g 9.18 9.18 9.69 8.62 9.31 N 2 1 R σ pπ 11.69 13.05 13.14 13.28 12.90 CO HCl HBr c Π : 3 u 1 A Π E 1 Π 1 A Π C 1 Π 1 A Π C 1 Π g * V : σ π R: σ 3pπ * C : π σ R : π σ * C : π σ R : π σ u 8.20 8.19 8.54 7.36 8.51 10.44 12.01 11.83 12.20 11.53 6.78 7.74 7.84 7.70 8.0 8.35 9.65 9.67 9.76 9.61 6.09 6.64 7.08 6.58 7.01* 7.30 8.41 9.29 8.57 8.74* D.P. Chong (Mol. Phys. 103, 749, 2005): 6.53 and 8.64 ev for A and C state of HBr by using SAOP (asymptotic V xc ) and augmented large basis set. * T e values used.

Next Implementation of MLR in Fourier-domain TDDFT

In Casida equation, 2 I ΩF = ω F I I Ω = δ δ ( ε ε ) ij, kl i, k j, l l k 2 + 2 ( f f )( ε ε ) K ( f f )( ε ε ) i j j i ij, kl k l l k If ( f f )( ε ε ) < 0, or ( f f )( ε ε ) < 0, what happens? i j j i k l l k ERROR: NaN 中 間 励 起 状 態 では 上 の 軌 道 の 占 有 数 は 下 の 軌 道 の 占 有 数 より 大 きいという 可 能 性 がある

Searching for Solution Non symmetric form of the Casida equation f f > 0 l k > 0 τ τ kτ lτ ε ε ε lτ εkτ > 0 εk τ ε lτ δστ, δik, δ jl, 2 Kij σ, klτ( ω) klτ fk τ fl τ 2 δστ, δik, δ jl, appl ω (Re δpkl τ )( ω) = δυij σ ( ω) ( fk τ fl τ)( εk τ εl τ)

lτ Further expressed as ε ε > 0 kτ klτ 2 δστ, δik, δ jl, ( εl τ εk τ) 2( fi σ f jσ)( ε jσ εi σ) Kij σ, klτ( ω) 1 + ω δ δ δ (Re δp )( ω) = δυ ( ω) 2 appl στ, ik, jl, klτ ijσ ( fi σ f jσ)( ε jσ εi σ) From the condition of poles, it satisfies that Ω Y Y 2 ( ω) I = ϖi I 2 Ω ijσ, klτ = δσ, τδi, kδ j, l ( εl τ εk τ) + 2( f f )( ε ε ) K Ω: 対 称 な 行 列 非 対 称 な 行 列 iσ jσ jσ iσ ijσ, klτ 数 学 的 には 厳 密 である This matrix is NOT symmetric (more time to be computed), but it should have the same eigenvalues as from the Casida equation. Moreover, the occupation number difference can be negative.

Singlet singlet transitions Straightforward calculation by the Casida equation in the spin unpolarized condition. ΩF = ϖ F 2 I I I ground intermediate excited 2 Ω ij, kl = δi, kδ j, l ( ε j εi ) + 2( f f )( ε ε ) K i j j i ij, kl

Implementation of MLR in the Casida Formalism Code: ABINIT Downloadable from http://www.abinit.org/ Open source GNU software project Plane wave Pseudopotential Hartwigsen Goedecker Hutter LDA & LSDA (Perdew Wang 92 functional)

ω I (ev) 16 Singlet-singlet 15 transitions 14 13 12 11 10 9 S4 8 S6 7 S8 6 0.0 0.2 0.4 0.6 0.8 1.0 q S3 S2 S5 S7 S1 Dependence of calculated vertical excitation energies on the fraction q of excited electron S1, S4, S8 (valence excitations): little dependence on q S2, S3, S5, S7 (Rydberg excitations): rapidly increase as q increases S6 (charge transfer excitation): rapidly increases as q increases At q=0 (ordinary LR), the Rydberg and charge transfer excitation energies are much underestimated.

Vertical excitation energies (in ev) Good Performance! State Transition Slater M-LR CI Exp. No. 1 a Π g V : σ g π g 8.57 9.20 9.69 9.31 S1 1 + N 2 a Σ R: σ g 3sσ g 12.25 12.34 12.17 12.20 S2 CO HCl g 1 c Π u R: σ g 3 1 A Π B 1 Σ + 1 A Π C 1 Π * V : σ π R: σ 3s * C : π σ R : π 4s pπ u 13.16 13.16 13.14 12.90 S3 7.32 8.22 8.54 8.51 S4 11.08 11.17 11.13 10.78 S5 7.69 7.87 7.84 8.0 S6 9.59 9.66 9.67 9.61 S7 1 C + * 5 Σu V : π π 3.01 6.35 6.35 5.88 S8 Slater: Slater transition state method ( SCF) M LR: midpoint value of modified LR

Summary (I) The modified LR works! Good performance on prediction of excitation energies of either valence, Rydberg or chargetransfer excitations. C. Hu, O. Sugino and Y. Miyamoto, Phys. Rev. A., 74, 032508 (2007). C. Hu and O. Sugino, J. Chem. Phys., 126, 074112 (2007).

Part IIの(2) TDDFT 線 形 応 答 理 論 に 基 づく 非 断 熱 結 合 係 数 の 計 算 Nonadiabatic couplings (NAC) from time dependent density functional theory within linear response

Photochemical Reactions Two central quantities: Potential energy surfaces Nonadiabatic couplings (NACs) Schematic photochemical reaction path (via conical intersection) to three different final structures: two products P and P and the starting reactant R Theor. Chem. Acc. (2006) 116: 87 105

Nonadiabatic quantum simulation Complete molecular wavefunction: (Nuclear wavefunction)(electronic wavefunction) 2 2 2 2 Φ(; rrt, ) = χn( Rt, ) Ψn(; rr), H = + V( R, r) 2 2 n 2M R 2m r 2 2 For fixed R: HR = + V( R, r), 2 HRΨ n(; r R) = En( R) Ψn(; r R) 2m r From the Schrödinger equation, we can get 2 2 i χm( R, t) = + E ( ) (, ) 2 m R χm R t t 2M R Nonadiabatic picture 2 M 2 2M n n Ψ m Ψ Born-Oppenheimer 近 似 : m R 2 R Ψ 2 n Ψ χn( R, t) R n χ ( R, t) n As the crossing point is approached, NACs Infinity! Nonadiabatic couplings (NACs) Hˆ Ψm Ψ n Ψ R m Ψ n = R E E 2 2 i χm( R, t) = + E ( ) (, ) 2 m R χm R t t 2M R n m 破 綻!

Nonadiabatic Transitions: Ubiquitous Conical intersections are ubiquitous in chemical systems Symmetry driven More often accidental Renner Teller glancing intersection Ideal conical intersections ( 円 錐 交 差 )(from wikipedia) Degenerate states at linear geometries Yarkony, Rev. Mod. Phys. 68, 985 (1996). Near crossings (or avoided crossings), nonadiabatic quantum simulation is necessary!

Requisites for nonadiabatic quantum simulation (1) Potential energy surfaces: E m excitation energies + total energy of the ground state (2) NACs near the crossing Ψ m R Ψ n, Ψ m 2 R 2 Ψ n When (1) and (2) are prepared, we can straightforwardly carry out nonadiabatic simulation (e.g. nuclear wavepacket propagation).

NACs from TDDFT? A ij, μ ( R) = Ψ i Ψ j R Many body wavefunctions are not available in DFT/TDDFT. For this reason, the computation of NACs has been regarded conceptually difficult for TDDFT. μ A real time TDDFT approach to compute NACs was presented by R. Baer, but it is also time consuming for small or medium sized molecular systems. R. Baer, Chem. Phys. Lett. 364, 75 (2002). New approach by TDDFT LR theory (in frequency domain) is necessary!

NACs from TDDFT: A rigorous formulation A ( R) I, μ 0 Hˆ Hˆ Ψ0 R Ψ μ I Ψ0 R Ψ μ = Ψ Ψ I = = R E E ω ˆ Using the perturbation operator ˆ H hμ R δv () t = hˆ ε () t μ appl μ μ I μ 0 I, we can introduce a perturbation as The dynamic polarizability from sum-over-states representation: α xz ( ω) = I 2ω Ψ hˆ Ψ Ψ hˆ Ψ I 0 x I I z 0 2 2 ωi ω I Chernyak and Mukamel, J. Chem. Phys. 112, 3572 (2000). The dynamic polarizability from frequency-domain TDDFT (Casida formalism): α xz + 1/2 + 1/2 2 hs FFS h Difference between ours and Casida formalism is only in the 2 2 perturbation operator: I ωi ω x I I z ( ω) = Hˆ R μ vs ˆR μ

NACs: S + 1/2 ˆ μ 0 hμ I 1/2 ωi Ψ Ψ = hs F A + 1/2 ( R) μ I, μ ω 3/2 I ijσ, klτ = = hs F Calculated by solving the Casida equation: δ δ δ στ, ik, jl, ( f f )( ε ε ) kτ lτ lτ kτ I I 規 格 化 条 件 : + I FF I = 1 ΩF = ϖ 2 F I I I Ω = δ δ δ ( ε ε ) ijσ, klτ σ, τ i, k j, l iσ jσ 2 + 2 ( f f )( ε ε ) K ( f f )( ε ε ) iσ jσ jσ iσ ijσ, klτ kτ lτ lτ kτ

To improve performance of LDA: MLR A + 1/2 μ I, μ ( R) ω 3/2 I = hs F I Nonadiabatic strength term Excitation energy The excitation energy is calculated from the response of the midexcited state, in which a half electron is promoted from the donor orbital to the acceptor orbital. The nonadiabatic strength term is not directly calculated from the mid-excited state, but from the response of the pure groundstate component of the mid-excited state. C. Hu, H. Hirai and O. Sugino, J. Chem. Phys. 127, 064103 (2007).

How to calculate the h matrix? The h matrix element is the recast of the explicit derivative of the many-body Hamiltonian with respect to nuclear coordinates in the Kohn-Sham density matrix. h Hˆ V = ϕ ϕ = ϕ ϕ en ij σ, μ i j i j R σ R σ μ μ Within the pseudopotential approximation: V = V = V + p E p l en pp loc l KB l l { } l l Vloc pl pl h, = R + R EKB p p E μ + μ KB Rμ ϕ ϕ ϕ ϕ ϕ ϕ ijσμ iσ jσ iσ l jσ iσ l jσ l In plane wave basis, the h matrix is calculated in the reciprocal space.

Implementation Code: ABINIT, a plane-wave pseudopotential program package (URL http://www.abinit.org) Exchange-correlation: LSDA Perdew-Wang 92 functional Ab initio norm-conserving pseudopotentials: Troullier-Martins (TM) Hartwigsen-Goedecker-Hutter (HGH) Fritz-Haber-Institute (FHI) pseudopotentials in TM scheme

A prototype system: H 3 Topic 1. Approaching conical intersection of H 3 : θ 0 in hyperspherical coordinates z 3 θ = 0 Equilateral triangle O 1 2 x 交 差

NAC (bohr -1 ) 交 差 点 に 近 づくと NAC が 発 散 する ρ = 2.5 bohrs, φ=120, θ = 10 0.1 200 150 100 50 x component of NACs Good agreement: MLR vs. Ref. (Abrol et al.) atom2_x,lr atom2_x,mlr Ref. NAC(bohr -1 ) JCP 115, 4640 (2001). NAC (bohr -1 ) 100 80 60 40 20 0 100 80 60 40 20 atom1_x, LR atom1_x, MLR Ref. 0 2 4 6 8 10 θ (degree) atom3_x,lr atom3_x,mlr Ref. 0 0 0 2 4 6 8 10 θ (degree) 0 2 4 6 8 10 θ (degree)

H Topic 2. Quantization of angular NACs in H 3 : ϕ q y x The Jahn-Teller effect H dx dy A12, ϕ = A12, x + A12, y dϕ dϕ = qa ( sinϕ + A cos ϕ) 12, x 12, y For the Jahn Teller ci, either perturbation theory or ab initio treatment, leads to a quantized (angular) NAC, namely A ϕ = 12, M. Desouter Lecomete, J. Phys. Chem. 89, 214 (1985). 1-1 (rad ) The result of ½ is limited to an infinitely small region surrounding the ci point. Usually the NACs will oscillate around ½. When the radius q is enlarged, the oscillation amplitude becomes larger. 2

Quantization of NACs in H 3 : Case 1 ϕ q 1.0 0.8 0.6 0.4 H H R HH =1.044 Å Radius q = 0.01 Å NAC (bohr -1 ) 0.2 0.0-0.2-0.4-0.6-0.8-1.0 0 60 120 180 240 300 360 ϕ (degree) LR MLR The MLR result is close to the quantized value (0.5), while that of LR is unreliable since the region is very close to the ci point.

Quantization of NACs in H 3 : Case 2 H ϕ q R HH =0.74 Å Radius q = 0.1 Å ϕ H NAC (rad -1 ) 1.0 0.8 0.6 0.4 0.2 0.0-0.2-0.4-0.6-0.8-1.0 0 60 120 180 240 300 360 ϕ (degree) LR MLR Ref. Oscillation around the quantized value becomes stronger; LR show clear improvement over LR results. Ref.: Halasz et al., J. Chem. Phys. 118, 3052 (2003).

H ϕ q R HH =1.044 Å Radius q = 0.3 Å Quantization of NACs in H 3 : Case 3 H NAC (rad -1 ) 1.0 0.8 0.6 0.4 0.2 0.0-0.2-0.4-0.6-0.8 LR MLR Ref. -1.0 0 60 120 180 240 300 360 ϕ (degree) Both LR and MLR give good results, with large amplitude of oscillation around the quantized value. By integrating NACs over the whole contour, the result is again quantized (π): Berry (geometric) phase effect can be studied. Ref.: Halasz et al., Chem. Phys. Lett. 358, 163 (2002).

2 1 A 12,ϕ (rad -1 ) Jahn-Teller 円 錐 交 差 in typical metal trimers ϕ q 1.4 1.2 1.0 0.8 y x 3 Ideal NACs as q = 0.02 bohr atom1-12.50-21.65 0.00 atom2 25.00 0.00 0.00 atom3-12.50 21.65 0.00 Quantization behavior reproduced =< Ψ Ψ > ϕ A12, ϕ 1 2 Theoretical value: 0.5 0.6 0.4 0.2 0.0 0 60 120 180 240 300 360 ϕ (degree) Angular NACs in Cu 3 trimer x y z Li 3 atom2 25.52 0.00 0.00 atom1-14.60-23.11 0.00 atom3-14.60 23.11 0.00 Na 3 atom2 24.25 0.00 0.00 atom1-13.42-21.28 0.00 atom3-13.42 21.28 0.00 K 3 atom2 24.53 0.00 0.00 atom1-13.38-21.51 0.00 atom3-13.38 21.51 0.00 Cu 3 atom2 25.13 0.00 0.00 atom1-13.97-22.38 0.00 atom3-13.97 22.38 0.00 Ag 3 atom2 26.67 0.00 0.00 atom1-14.13-22.88 0.00 atom3-14.13 22.88 0.00 Au 3 atom2 26.62 0.00 0.00 atom1-13.58-22.32 0.00 atom3-13.58 22.32 0.00

Renner-Teller intersections x, y and z components of NACs when q = 0.1 bohr BH 2 atom1-5.43 0.00 0.00 atom2 12.94 0.00 0.00 atom3-5.43 0.00 0.00 AlH 2 atom1-6.29 0.00 0.00 atom2 13.66 0.00 0.00 atom3-6.29 0.00 0.00 CH + 2 atom1-5.82 0.00 0.00 atom2 13.12 0.00 0.00 atom3-5.82 0.00 0.00 SiH + 2 atom1-6.51 0.00 0.00 atom2 10.51 0.00 0.00 atom3-6.51 0.00 0.00 H 1 r y 2 q Ideal NACs as q = 0.1 bohr atom1-5.00 0.00 0.00 atom2 10.00 0.00 0.00 atom3-5.00 0.00 0.00 Reasonable agreement with theoretical model x r H 3

Challenges of the pseudopotential approximation in the computation of NACs The x components of NACs for NH 2 and H 2 O + molecules NH (a) 2 NH (b) 2 H 2 O + (a) H 2 O + (b) Model H1-5.09-5.31-5.22-5.33-5.00 X 41.28 13.20 41.34 12.58 10.00 H2-5.09-5.31-5.22-5.33-5.00 (a) Nonlocal (b) local pseudopotentials A h S FI ( ), + 1/2 μ I, μ R = 3/2 ωi ijσ, μ Hˆ V ϕ ϕ = ϕ ϕ en iσ jσ iσ jσ Rμ Rμ Enormous freedom in constructing pseudopotentials h Difference of H and X atoms only comes from V en (V PP ). = C. Hu, H. Hirai and O. Sugino, J. Chem. Phys. 128, 154111 (2008).

All electron calculation of NAC by TDDFT Difference of basis sets Pseudopotential / Plane wave: h H V PP ele-nucl mnσ, μ = ϕm σ ϕn σ = ϕm σ ϕn σ = ϕm σ ϕn σ Rμ Rμ Rμ All electron / atomic orbital basis: H H h = ϕ ϕ = C χ C χ mnσ, μ mσ nσ imσ i jnσ j Rμ i Rμ j = i j C C χ H χ imσ jnσ i j Rμ Effects of atomic orbital basis must be considered in the all electron scheme when constructing h V

Molecular Hamiltonian in atomic orbital representation using second quantization 1 H = t c c + ij kl c c cc ij ij i j i j l k 2 ijkl r r r r () ZI v nucl r = r R * 1 2 tij = χi () r + vnucl() χ j () d 2 * * χi ( r1) χk( r1) χ j( r2) χl( r2) ij kl ( ik jl) = drdr r r 1 2 1 Tommasini, Chernyak and Mukamel, Int. J. Quant. Chem. 85, 225 (2001). 2 Szabo and Ostlund, Modern quantum chemistry, 1989. I 1 2 3 Orthonormal basis set is used in the expression Tommasini et al. derived the expression of h: R H x x x h = h ij c i c j I ij h t V ρ x ij x = + ij kl ee x kl I

h matrix element Hartree Fock (HF) exact exchange only: h = t + ( ij kl) ρ ( ik jl) x, α x x tot x α ij ij kl kl kl kl General exchange correlation (xc) functionals: h = t + ( ij kl) ρ x, α x x tot ij ij kl kl x αα α x αβ β + ( ij) fxc kl ρkl + ( ij) fxc kl ρkl kl kl αα x α αβ x β + ij fxc ( kl) ρkl + ij fxc ( kl) ρkl kl kl xc kernel: f αβ xc 1 (, rr ) = ρ 2 I tij = i j 2 I R I 2 δ Exc[ ρα, ρβ ] δρ () r δρ ( r ) α r β Z

Orthonormal Transformation From atomic orbital basis χ, usually nonorthonormal, to orthonormal basis χ : χ = χs 1/2 The molecular orbital expansion coefficients are transformed as: 1/2 C= S C The density matrix is transformed as: ρ = S ρ S 1/2 1/2 ij ik kl lj kl The matrix of one-electron operator is: S: Overlap matrix x x t = S ts + S t S + S t S ( 1/2 ) 1/2 1/2 x 1/2 1/2 ( 1/2 ) S = χ χ ij i j x

Our method for derivative of S -1/2 SV = Vε, S = Vε V, S = Vε V 1/2 1/2 1/2 1/2 1/2 1/2 1/2 1/2 ( R ) + ( R ) = R S S S S S 1/2 1/2 1/2 1/2 ( R ) + ( R ) = R S Vε V Vε V S S V ( S ) Vε + ε V ( S ) V = V ( S) V 1/2 1/2 1/2 1/2 R R R 1/2 1/2 1/2 1/2 [ V ( RS ) V] ij[ ε ] jj + [ ε ] ii[ V ( RS ) V] ij = [ V ( RS) V] ij 1/2 [ V ( RS ) V] ij = [ V ( RS) V] ij 1/2 1/2 [ ε ] ii + [ ε ] jj S = V{[ V ( S ) V] } V 1/2 1/2 R R ij S = S ( S ) S R 1/2 1/2 1/2 1/2 R

Implementation Code: GAMESS (Open source quantum chemistry package) http://www.msg.chem.iastate.edu/gamess/ TDOSCALC, UTDOSCALC NAC: Analogy to oscillator strength h:analogy to force Minimization of coding by modifying existing subroutines HD_HELFEY HD_TVDER HD_JKDER HD_DFTDER Hellmann Feynman part: ( r R) μ χi Z χ 3 j r R Kinetic and V nuclear contribution via atomic orbital basis : χ χ i j χ + χ + χ + χ R T V T V j nucl j i nucl j μ Rμ Contribution of the two electron operator via atomic orbital basis: ρ χ χ 1 χ χ kl i j k l r r kl R μ +... Contribution of the xc part via atomic orbital basis ( 現 段 階 HFのみ)

Calculation Results Renner Teller systems: q=1.0 bohr Tab1 x component: Reasonable agreement H 1 r y X q x r H 3 Strict Satisfaction of sum rule in TD/HF Pulay terms

Jahn Teller systems Hyperspherical coordinates: ( ρ, θ, ϕ) 3 H 3 NAC (bohr -1 ) 24 20 16 12 8 4 TD/HF TD/LDA Ref. 1 1 2 LDA gives better performance than HF when θ is small. 2 y ϕ q x 3 q = 0.5 bohr Li 3 : exchange correlation effects (more than HF exchange) might be important 0 2 4 6 8 10 θ (degree) TD/HF CASSCF H 3 Atom 1 0.468 0.467 Atom 2 0.937 0.918 Atom 3 0.468 0.467 Li 3 Atom 1 0.158 0.487 Atom 2 0.315 0.975 Atom 3 0.158 0.487 Hu, Sugino, Tateyama, J. Chem. Phys. 131, 114101 (2009).

Summary (II) 1. 多 体 波 動 関 数 で 定 義 した 非 断 熱 結 合 係 数 (NAC)に 対 して TDDFT 線 形 応 答 理 論 に 基 づき 多 体 波 動 関 数 を 求 めずNACの 厳 密 な 計 算 式 を 導 出 した 2. 平 面 波 擬 ポテンシャル 枠 組 みで 実 証 したところ 様 々の 系 で 計 算 精 度 が 良 好 であることを 確 認 した 一 方 擬 ポテンシャル 近 似 の 問 題 点 も 明 らかになった 3. さらに 全 電 子 原 子 基 底 計 算 手 法 を 開 発 し TDDFTによるNACの 計 算 の 妥 当 性 を 調 べた その 結 果 Hartree Fock exact exchangeを 用 いた 計 算 結 果 は NACのsum ruleをよく 満 たすことが 分 かった Very recent work: Accurate calculation of NAC in the pseudopotential planewave TDDFT (to be submitted) Work on the way: 2 nd order NAC from TDDFT ( 日 本 物 理 学 会 秋 季 大 会 で 発 表 する 予 定 )

謝 辞 杉 野 修 准 教 授 ( 東 大 物 性 研 ) 宮 本 良 之 博 士 (NEC) 館 山 佳 尚 博 士 ( 物 質 材 料 研 究 機 構 ) 平 井 宏 俊 さん( 東 大 物 性 研 杉 野 研, D3) 渡 辺 一 之 教 授 ( 理 科 大 )