10_excited - YoshioNishimoto/sandbox GitHub Wiki

励起状態計算

一言で励起状態計算と言ってもいろいろな方法がありますが、 ここでは単一スレーター行列式（単配置・単参照）を用いた計算手法により励起状態のエネルギーを計算する二つの手法を取り扱います。 具体的には、(linear-response) time-dependent density functional theory (TD-DFT)と、 equation-of-motion coupled-cluster (EOM-CC)を使ってみます。 どちらの手法も単配置・単参照ということで、特に面倒な入力ファイルの作成が必要ありません。

ここで扱う励起状態は、回転や振動の励起状態ではなく、電子の遷移に伴う励起状態です。 すでに授業でも勉強したと思いますが、例えばHOMOの電子がLUMOへ遷移するような過程です。 多くの場合、励起エネルギーは数eV単位となり、可視光領域の光を吸収します。 これらの計算で吸収スペクトルのシミュレーションや、蛍光・燐光のシミュレーションなどを行う事ができます。

TD-DFT計算

TD-DFT計算は、configuration interaction singles (CIS)計算と類似の点が多いです。 とは言え、最近はTD-DFTはCISの計算コストとほとんど同じで、なおかつTD-DFTの方が精度が良いため、 CISが用いられることは少ないです（個人の感想です）。 なので、とりあえずはTD-DFTが扱えれば励起状態計算がある程度できると主張しても良いと思います。 TD-DFTとCISの類似性に関しては、TD-DFTからDFTの代わりにHartree--Fockを用いて、さらにTamm--Dancoff近似を用いると、CISになります。

とりあえず、以下の通りethyleneの励起状態計算をしてみます。 最初に一重項基底状態（S0と言います）で構造最適化して、励起状態計算を行うというのが良くある流れです。 これで、分子が吸収するエネルギー（垂直励起エネルギーと呼ばれることが多い）や光の波長を計算することができます。 以下の構造は既にS0で構造最適化を行ってあります。

  1 from pyscf import gto, dft, tddft
  2 
  3 mol = gto.M(verbose=4,
  4             atom= [["C",(0.00000000, 0.00000000,-1.25818248)],
  5                    ["C",(0.00000000, 0.00000000, 1.25818248)],
  6                    ["H",(0.00000000, 1.74600192,-2.34319393)],
  7                    ["H",(0.00000000,-1.74600192,-2.34319393)],
  8                    ["H",(0.00000000, 1.74600192, 2.34319393)],
  9                    ["H",(0.00000000,-1.74600192, 2.34319393)]],
 10             basis="6-31G(d)",unit="Bohr")
 11 mf = dft.RKS(mol)
 12 mf.xc = 'b3lyp'
 13 mf.kernel()
 14 
 15 mytd = tddft.TDDFT(mf)
 16 mytd.kernel()
 17 mytd.analyze()

11から13行目で基底状態でのSCF計算を行います（必須）。 TD-DFT計算は、基底状態の電子構造（分子軌道係数）を用いて線形近似により励起エネルギーを計算するからです。 15行目でTD-DFT計算を行う準備をして、16行目で実際の計算を行います。 そして17行目で簡単な解析を行います。 つまり、新しく最後の三行（とimport）を加えるだけでTD-DFT計算ができます。

SCF計算の箇所は置いておいて、TD-DFT計算の結果以下の出力を得られます。

 99 ******** <class 'pyscf.tdscf.rks.TDDFT'> for <class 'pyscf.dft.rks.RKS'> ********
100 nstates = 3 singlet
101 wfnsym = None
102 conv_tol = 1e-09
103 eigh lindep = 1e-12
104 eigh level_shift = 0
105 eigh max_space = 50
106 eigh max_cycle = 100
107 chkfile = /data/home/nisimoto/lib/pyscf/test/test/excited/tmpnfe0x3qq
108 max_memory 4000 MB (current use 101 MB)
109 
110 
111 Excited State energies (eV)
112 [8.30427957 8.4707365  9.23016164]
113 
114 ** Singlet excitation energies and oscillator strengths **
115 Excited State   1:      8.30428 eV    149.30 nm  f=0.3512
116        8 -> 9        -0.70408
117 Excited State   2:      8.47074 eV    146.37 nm  f=0.0000
118        7 -> 9         0.70162
119 Excited State   3:      9.23016 eV    134.33 nm  f=0.0003
120        8 -> 10        0.70537
121 
122 ** Transition electric dipole moments (AU) **
123 state          X           Y           Z        Dip. S.      Osc.
124   1        -0.0000     -0.0000     -1.3139      1.7264      0.3512
125   2         0.0000      0.0000      0.0000      0.0000      0.0000
126   3        -0.0333     -0.0000      0.0000      0.0011      0.0003
127 
128 ** Transition velocity dipole moments (imaginary part, AU) **
129 state          X           Y           Z        Dip. S.      Osc.
130   1        -0.0000      0.0000     -0.4075      0.1661      0.3628
131   2         0.0000     -0.0000     -0.0000      0.0000      0.0000
132   3        -0.0356      0.0000     -0.0000      0.0013      0.0025
133 
134 ** Transition magnetic dipole moments (imaginary part, AU) **
135 state          X           Y           Z
136   1         0.0000     -0.0000      0.0000
137   2        -0.0000      0.0000      1.3958
138   3         0.0000      0.0000      0.0000

一番重要なのは、mytd.analyze()により得られる114から120行目でしょう。 励起エネルギー、そのエネルギーに対応する波長、振動子強度（f; oscillator strength）、 （おそらく）その励起状態で主要な電子配置が示されています。 例えば第一励起状態（Excited State 1:; S1）に関して言えば、 励起エネルギーが8.304 eV、波長が149.30 nm(\approx 1240/8.304)、振動子強度が0.3512、 最も主要な電子配置が8番目の軌道から9番目の軌道へと電子を励起する電子配置であると書いてあります（これらはどのような軌道でしょうか？）。 実験によると160から165 nm辺りの吸収があるため、やや励起エネルギーを過大評価している感じがします（が、こんなもんです）。

振動子強度は遷移のしやすさを示す指標です。 実験で得られる吸収スペクトルには様々なピークが現れます。 この高さに対応するのが振動子強度です。 f=0.3512は他の状態に比べれば割と大きく、この遷移が吸収スペクトルで高い吸収強度として現れそうだということが言えます。 実際に実験結果と比較する場合には、ガウス関数やローレンツ関数を用いて広がりのあるピークとして表現して、近接する励起状態と重ね合わせたりすることが多いです。

TD-DFT計算の詳細に関しては、出力ファイルには何も書いてありませんが、大きな行列を対角化しています（Casida's equation）。 ethyleneのような小さな分子であればともかく、大きな分子では対角化すべき行列の要素（普通のTD-DFTだと、確か4Nocc2Nvirt2の次元）を全てメモリに格納することができません。 そこで、Davidson Methodというアルゴリズムを用います。 このアルゴリズムは、大きな行列をそのまま対角化するのではなく、 すべての固有値のうち固有値の小さいいくつか（デフォルトでは3状態ですが、例えばmytd.nstates=5を入れておくと5状態を計算します）の固有値・固有ベクトルのみを計算します。 行列要素をメモリに格納することなく、大きな行列とベクトルのかけ算（ただし行列は保存せず、on-the-flyでかけ算を行います）・小さな部分空間（Krylov subspace）での対角化のみで計算できるので、大きな分子での励起状態計算が可能になります。

EOM-CC計算

同じ分子を使ってEOM-CCの計算をしてみます。 より具体的にはEOM-CCSDと呼ばれる手法で、単参照電子相関理論で出てきた結合クラスター理論を用いて一電子・二電子励起を考慮しながら励起状態を取り扱っていると思われます。 なので、精度的にもCCSDレベルとなります。

  1 from pyscf import gto, scf, cc
  2 
  3 mol = gto.M(verbose=4,
  4             atom= [["C",(0.00000000, 0.00000000,-1.25818248)],
  5                    ["C",(0.00000000, 0.00000000, 1.25818248)],
  6                    ["H",(0.00000000, 1.74600192,-2.34319393)],
  7                    ["H",(0.00000000,-1.74600192,-2.34319393)],
  8                    ["H",(0.00000000, 1.74600192, 2.34319393)],
  9                    ["H",(0.00000000,-1.74600192, 2.34319393)]],
 10             basis="6-31G(d)",unit="Bohr")
 11 mf = scf.RHF(mol)
 12 mf.kernel()
 13 
 14 mycc = cc.RCCSD(mf)
 15 mycc.kernel()
 16 e_ip, c_ip = mycc.ipccsd(nroots=3)
 17 e_ea, c_ea = mycc.eaccsd(nroots=3)
 18 e_ee, c_ee = mycc.eeccsd(nroots=3)

14・15行目で基底状態のCCSD計算を行い、 16・17・18行目でionization potential (IP; 電子を一つ取り除く)、electron affinity (EA; 電子を一つ追加する)、electron excitation (EE; 電子励起をする)のEOM-CCSD計算を行います。 ここでは、最後の行のe_ee, c_ee = mycc.eeccsd(nroots=3)の結果が、上のTD-DFT計算のように励起エネルギーの計算を行う行に対応します（が、tripletやspin-flipの結果が含まれるので、よく分かりません）。

とりあえず確認するべきは、一重項励起のエネルギーでしょう。

143 EOM-CCSD root 0 E = 0.3391302004691231  qpwt = 0.974217  conv = True

単位はHartreeなので、9.228 eVぐらいになります。 こちらも割と過大評価するようです。

EOM-CCSDは、日本の某研究室で開発されたsymmetry-adapted cluster CI (SAC-CI)という手法と等価であると言われています。 どちらの名前を使うかは何とも言えませんが、用いたプログラムで使われている名前を書くのが一般的です。

基本的にEOM-CCの計算はTD-DFTに比して計算コストがかなり高いです。 その上、TD-DFTでもかなり良い精度で励起状態計算ができるため、特に最近は励起状態計算と言えばTD-DFTというぐらいです。 しかし、TD-DFTは汎関数依存性があるため、場合によっては大きく実験から外れる場合があります。 有名なのは、電子移動が伴う励起（例えばdonor-acceptor complexでの励起）は、長距離補正を行った汎関数を用いる必要があります。 ですが一方で、長距離補正を行った汎関数は、π-π*励起のエネルギーをかなり過大評価しやすいです。

励起状態で構造最適化をすることも可能です。 詳細は省きますが、例えばB3LYP汎関数を用いたTD-DFT計算（TD-B3LYPと書かれることが多い）では

mf = dft.RKS(mol)
mf.xc = 'b3lyp'
mytd = tddft.TDDFT(mf)
mol_eq = optimize(mytd)

でできるようです。 励起状態で構造最適化し、その極小値での垂直励起エネルギーは、（一重項の場合）蛍光のエネルギーに相当します。 三重項で同じ計算を行えば、燐光に相当します。

上記計算はすべて真空中での計算です。 溶媒効果を取り入れてTD-DFT計算を行うことも可能（PySCFができるかは分かりませんが）ですが、正しく自由エネルギーを計算するのはなかなか複雑です。

演習問題

  • ブリルアンの定理によると、「（HFあるいはDFTといった）1粒子法によって解かれた基底電子状態がすでに1電子励起配置と基底状態配置の配置間相互作用を暗黙的に含む」とされています。ですが、TD-DFTやCISは一電子励起配置のみを扱う手法にもかかわらず、（当然ながら）励起エネルギーはゼロではない値となります。基底状態と励起状態で、一電子励起配置がどのようにエネルギーに寄与するかを（できればブラケット表記で）説明してください。Szabo: Section 4.1, 4.2あたりが基底状態の参考に、Q-Chemのページが励起状態の参考になると思います。
  • 適当な分子を使って、励起エネルギーにどの程度の汎関数依存性が出るかを検証してみましょう。
  • 上の分子をcc-pVTZ基底関数で計算して、実験で得られる吸収スペクトルの波長（本文参照）と比較してみましょう。 EOM-CCSD/cc-pVTZはそれなりに時間がかかります。
  • TD-B3LYP/6-31G(d)を用いて、ethyleneのS1状態から輻射緩和されるときに発する光のエネルギー（または波長）を計算してみましょう。また、断熱励起エネルギーを計算してみましょう。断熱励起エネルギーは、基底状態の安定構造のエネルギーと、励起状態の安定構造のエネルギー差（ES1 - ES0）で定義します（ちなみにゼロ点振動補正は必要ないです。ゼロ点振動補正を加えたエネルギー差は、0-0励起エネルギーと呼ばれるものに対応します）。
  • TD-B3LYP/6-31G(d)を用いて、できるだけ垂直励起エネルギー（S0状態での構造最適化を含む）の値が小さくなる分子を設計しましょう。第3周期以降の元素を用いる場合は、計算結果を保証するデータや文献を提示できるようにしましょう。
  • 任意の汎関数や基底関数を用いて、TD-DFTでtrans/cis-azobenzeneの垂直励起エネルギーや蛍光エネルギーを計算してください。 適当な文献で実験値または高精度計算手法の値を得て、参考文献を示しつつ値を比較してください。azobenzeneはC12H10N2ということでそれなりに分子が大きいので注意してください。
⚠️ **GitHub.com Fallback** ⚠️