虚数時間発展法による電子状態計算手法の開発

Abstract

We have been studying the time-evolution technique in imaginary time. The electronic structure calculated in this technique is valid for the ground state, which is easily confirmed from the completeness of the wavefunction. However, the excited states can be also obtained by Gram-Schmidt orthogonalization technique. As a practical 3D calculation, we applied this technique to hydrogen atom and molecule systems. The results are in good agreement with their analytical solutions, including the excited states. The purpose of this study is to present the effectiveness of this technique for practical electronic structure calculations. To realize it, we will start to apply this to hydrogen and practical small molecule systems.

Translated Abstract

We have been studying the time-evolution technique in imaginary time. The electronic structure calculated in this technique is valid for the ground state, which is easily confirmed from the completeness of the wavefunction. However, the excited states can be also obtained by Gram-Schmidt orthogonalization technique. As a practical 3D calculation, we applied this technique to hydrogen atom and molecule systems. The results are in good agreement with their analytical solutions, including the excited states. The purpose of this study is to present the effectiveness of this technique for practical electronic structure calculations. To realize it, we will start to apply this to hydrogen and practical small molecule systems.

1 序論

光学的な物性を探るためには電子の励起状態を求める必要がある．そのための手法として時間依存密度汎関数法(TDDFT)などが用いられる．特に実時間を用いる方法が数値計算上も安定に計算することができる．この実時間のTDDFTを用いる電子状態計算では，時間発展の前に，まず定常状態(基底状態)を求める必要がある．その状態に外部刺激を摂動として加え実時間でその応答を追跡するのが基本である．また十分収束させておかなければ，その後の時間発展の追跡に支障をきたす．また実時間のTDDFTでは占有軌道のみの電子の動的な変化から光学的な特性などを算出する [1]．

この初期の定常状態を求める際にも同じ時間発展の枠組で計算が可能である．基底状態から逐次の状態を求めることができる虚数時間発展法を使うことである．この方法は原子・原子核などの分野では古くから用いられている実績のある方法であるが，すべての計算を実時間でのTDDFTと同じ時間発展の枠組みの中で計算でき，実時間発展においてで開発されてきた同じ計算技術が使えるので非常に効率が良い．これを示すことが目的である．また，通常の固有状態計算法より，時間発展計算の方が並列化効率よく，全計算量が少なくて済むということが知られている [1, 2]．

解析手法の確認は固有状態の解析解が良く知られている系を用いた．この手法を用いて水素原子および水素分子の固有状態の解析を行った．

2 計算手法

多体電子系の電子状態を調べるため，時間に依存するKohn-Sham方程式において，実時間を虚数時間に置きなおす($\tau =-it$) と，虚数時間発展の方程式が得られる [3]．これは式(1)で表されるように拡散方程式の形をしており，数値的にも安定に計算できることが知られている．   

$-\frac{\partial }{\partial t}\psi \left(r,\tau +\Delta \tau \right)=\left[-\frac{1}{2}{\nabla }^2+V\left(r\right)\right]\psi \left(r,\tau \right)$(1)

ここで $V\left(r\right)=V_{ps}\left(r\right)+V_H\left(\rho \left(r\right)\right)+V_{xc}\left(\rho \left(r\right)\right)$である．$V_{ps}$は擬ポテンシャルであり，$V_H$はハートレー項，$V_{xc}$は交換相関相互作用ポテンシャルである．この解は一般的に式(2)で表すことができる．   

$\psi \left(r,t\right)=\sum _{k=1}^{\infty }{C}_k{\varphi }_k\left(r\right)\text{ exp}\left(-{\epsilon }_k\tau \right)$(2)

ここで ${\varphi }_k$ と ${\epsilon }_k$ は$k$ 番目の各固有状態と固有エネルギーを表わす．虚数時間$\tau$を発展させていくと，   

$\begin{array}{c}{e}^{-H\tau }\psi \left(r,\tau \right)\text{=}\sum _n{e}^{-E_n\tau }|{\varphi }_n⟩⟨{\varphi }_n|\psi ⟩\\ =e^{-E_0\tau }\sum _n{e}^{\left(E_0-E_n\right)\tau }|{\varphi }_n⟩\text{}\text{}\text{}⟨{\varphi }_n|\psi ⟩\\ \to e^{-E_0\tau }\left|{\varphi }_0⟩⟨{\varphi }_0\right|\psi =e^{-E_0\tau }|{\varphi }_0⟩\left(\tau \to \infty \right)\end{array}$(3)

固有エネルギーの差が指数関数的に減少するためエネルギーが一番低い基底状態が求まる．これは波動関数の完全性$1={\Sigma }_n|{\varphi }_n⟩⟨{\varphi }_n|$ から次のようにして説明できる．この式(3)により，虚数時間発展をすると，基底状態より高い固有エネルギーを持つ状態は指数関数的減衰するので，固有状態 ${\psi }_0\left({\epsilon }_0\right)$ しか残らない．十分な時間発展を行うと精度よく，確実に基底状態を求めることができる．   

$\underset{\tau \to \infty }{\lim }\psi \text{(}r,\tau \text{)}\approx C_0{\psi }_0\left(r\right)\text{exp}(-{\epsilon }_0\tau \text{)}$(4)

この解を一例としてCrank-Nicolson法を使って求めた．実時間発展でもよく用いられる方法である．時間を現時間$\tau$から $\nabla \tau /2$ だけ進めて，求める状態から $\nabla \tau /2$ だけ戻すと，そこで状態は等しいことから，   

${\text{e}}^{\frac{\Delta \tau }{2}H}\psi \text{(r},\tau \text{+}\Delta \tau \text{)}={\text{e}}^{\frac{-\Delta \tau }{2}H}\psi \text{(r},\tau \text{)}$(5)

となる．この方法によって，虚数時間発展の演算はテーラー展開の1次まででもユニタリー性を確保でき，計算の安定性を保てる.

もちろん虚数時間発展法では電子状態の基底状態しか求められない．他の固有状態を求めるには，他の状態との直交性から， Gram-Schmidtの直交化法を用いて求める [3]. 第$n-1$番目の状態までわかっているとすると，第$n$番目の状態は，次のように逐次的に求めることができる．   

$\begin{array}{c}\psi {\text{'}}_n={\tilde{\psi }}_n-\sum _{i=1}^{n-1}\frac{1}{{\left|{\psi }_i\right|}^2}⟨{\psi }_i|{\psi }_n>{\psi }_i\\ {\psi }_n\text{ = }N\text{(}\psi {\text{'}}_n)\end{array}$(6)

ここで$N$は規格化を表わす．

3 結果と考察

虚数時間発展法のプログラムでの検証は，1次元井戸型ポテンシャル，モースポテンシャルのなどの基本なポテンシャルおよび水素原子など解析解のある単純な系にて計算を試みた．いずれの計算においても解析解と同じ値を得られた．この手法を，対称性を用いない3次元に拡張しTDDFTで実時間発展を行うものと同じグリッド${10}^3 a.u.\left(\Delta x,\Delta y,\Delta z=0.1 a.u.\right)$ を用いて，水素原子および水素分子の解析を行った．

Figure 1は水素原子の固有状態を求めた結果である．もちろん水素原子の場合，通常の手法$\psi \left(r\right)=R\left(r\right)Y\left(\theta ,\phi \right)$とし，動径方向の波動関数を更に$rR\left(r\right)=F\left(r\right)$と変換して求める手法でも求めることができるが，いずれの場合も動作確認が目的であるため，単純な軌道角運動量 $l=0$ の場合のみを用い，解析解と一致していることを確認した [3, 4]．Figure 2は水素分子の結合状態と反結合状態を示したものである．基底状態である結合状態を虚数時間発展法により求め，励起状態である反結合状態は十分収束させた基底状態にGram-Schmidtの直交化を行い求めた．いずれの場合も様々な材料系での利用を考慮しTroulier-Martinsの擬ポテンシャルを用い，境界条件として遠方での漸近形はSlater型を採用した[5]．

Figure 1.

 The$$orbitals of hydrogen atom

Figure 2.

 The$$orbitals of hydrogen molecule of (a) bonding and (b) anti-bonding states

実時間のTDDFTを開始するには，ここで得られた占有状態に外部刺激を摂動として加えその応答を追跡する．もちろん，虚数時間発展による計算は，実時間発展による計算のための初期値を作る計算であり，TDDFT計算とは直接関係がない．しかし時間発展に関しては，虚数時間にせよ実時間にせよ，同じ枠組で同様の操作を行うので効率よく計算を移行することが可能である．

虚数時間発展法では，基底状態以外の状態が指数関数的に減衰するため効率よく基底状態を求めることができるが，今後，小型の分子など様々な系での比較を行う予定である．

4 まとめ

電子の固有状態の解析に虚数時間発展法を用いた．適切な遠方での漸近形を与えれば，水素原子および水素分子についてではあるが，効率よく確実に固有状態を求めることができる．また電子状態計算のプログラム全体を同じ時間発展の枠組みの中で計算できるという利点がある．

謝辞

この研究はJSPS科研費16K05047 の助成支援を受けたものである．

参考文献

• [1]Y. Zempo, N. Akino, M. Ishida, E. Tomiyama and H. Yamamoto, J. Phys.: Conf. Ser. 640, 0l2066 (2015).
• [2]   A. K. Roy, J. Math. Chem., 52, 2645 (2014). doi:10.1007/s10910-014-0407-0
• [3]   G. Kresse, J. Furthmüller, Phys. Rev. B, 54, 11169 (1996). doi:10.1103/PhysRevB.54.11169
• [4]   F. Hirono, M. Iwasawa, S. S. Kano, Y. Zempo, J. Comput. Chem. Jpn., 18, 159 (2019).
• [5]   K. Takayanagi, Chapter 2. "Genshi Bunshi Butsurigaku", Asakura Pulishing Co., Ltd, 2000.