4P122

Reptation Monte Carlo 法を用いた多原子分子の高精度振動状態解析(横浜市大院・生命ナノ)○藤岡 蔵、北幸海、立川仁典

Accurate vibrational state analysis of polyatomic moleculeswith Reptation Monte Carlo method(Yokohama City University) ○ Osamu Fujioka, Yukiumi Kita, Masanori Tachikawa【序論】 赤外分光法やラマン分光法は、分子の構造決定や分子内・分子間相互作用の解析に用いられており、現代の分子科学において必要不可欠な実験法である。例えば赤外分光法では、対象となる分子の赤外吸収スペクトルを測定し、そのスペクトルを 読み解く ことで分子の構造情報を抽出している。しかし一般的に、測定によって得られるスペクトルは極めて複雑であり、スペクトル情報を分子の構造情報へと還元するためには、計算化学的手法によるスペクトルの帰属が必要不可欠となる。例えば、分子軌道法に代表される第一原理計算を用いて、分子の平衡構造における基準振動解析から振動スペクトルを予測し、その結果を実験と比較することで、スペクトル情報と分子の構造情報を対応付ける。しかしながら、一般的に良く用いられる基準振動解析では、ポテンシャルの非調和性や振動モード間のカップリングを無視するため、振動数を過大評価してしまう問題が知られている [1−3] 。このような問題を解決するために、当研究室ではこれまでに、量子モンテ
カルロ (QMC) 法の一つである Reptation Monte Carlo (RMC) 法 [4] に、Multi-Product (MP)展開法 [5] を組み合わせた、新しい量子モンテカルロ法 (RMC+MP 法) の開発を行ってきた [6] 。
本研究では、ポテンシャルの非調和性や振動モード間のカップリングを精密に考慮した振動状態における物理量の精密算定および振動状態における分子構造の解明を目的に、RMC+MP 法による振動状態解析の精度検証を行った。具体的には、H2 O とその同位体である HOD、D2 O の分子振動状態へ RMC+MP 法を適用し、基本振動数、存在確率密度、回転定数の解析を行った。
【計算方法】 時間依存 Schrödinger 方程式の虚時間発展を利用した QMC 法では、以下のように、任意の試行波動関数 ΨT に虚時間プロパゲータを作用させることにより、系の正確な固有状態 Ψ0 を抽出する。
β|Ψ0 ⟩ = lim e− 2 Ĥ |ΨT ⟩(1)β→∞

ここで β は虚時間、Ĥ はハミルトニアンを表す。RMC 法では、仮想的な分配関数

Z0 ≡ ⟨Ψ0 |Ψ0 ⟩ = lim ⟨ΨT |e−βĤ |ΨT ⟩β→∞(2)を導入し、十分長い虚時間後 (β → ∞) の密度分布をサンプリングすることにより、原理上、系の正確な波動関数とエネルギー固有値を求めることができる。具体的には、式 (2)
中の虚時間プロパゲータを短時間分解した後、短時間近似によって得られる分配関数を用いる。例えば、RMC 法の二次の短時間近似を用いた分配関数は以下のように表される。
 N−1∫∫N−1∏∏dτ
Z0 =· · · ΨT (x0 ) g(xi , xi+1 ; dτ) ΨT (xN ) ×(3)e− 2 {ε(x j )+ε(x j+1 )} dx0 · · · dxN
i=0j=0ここで N は虚時間プロパゲータの分割数、dτ(= β/N) は虚時間ステップを表す。式 (3) のg は |ΨT |2 に従う分布を生成する短時間プロパゲータ、ε(x) ≡ Ψ−1T (x) ĤΨT (x) で表される局所エネルギーを含む項は、生成したパスの重み因子と解釈できる。
物理量の解析方法について、まず全エネルギーは以下の mixed estimator を用いて算出した。
⟨ΨT |Ĥe−βĤ |ΨT ⟩ ⟨ΨT |Ĥ|Ψ0 ⟩E0 = lim=(4)β→∞ ⟨Ψ |e−βĤ |Ψ ⟩⟨ΨT |Ψ0 ⟩TT
存在確率密度等のハミルトニアンと非可換の物理量は、虚時間パスの中央が正確な分布|Ψ0 |2 に対応することから、その配置を用いて算出した。
【計算の詳細】 基準振動座標と多次元ポテンシャルエネルギー曲面 (PES) は、分子軌道計算プログラム Gaussian09 を用いて、CCSD(T)/aug-cc-pVTZ レベルで算出した。配置空間中の grid 点の数は 9261(= 213 ) 点である。QMC 計算における試行波動関数には平均場近似に対応する以下の VSCF 形式を用いた。
Nbasis (( (i)))ν2∏ ∑ (i) (i) αn (Qi − Q(i)c )(i)(i) (i)Ψ(Q; C, α, Qc ) =cn φn  , φn = Nn Hn exp −(5)2i=1 n=0
ここで Nn(i) は規格化因子、Hn(i) はエルミート多項式を表す。本研究では展開係数 c(i)n 、軌道(i)(i)[7]
指数 αn 、基底関数中心 Qc を変分パラメータとし、線形最適化法 により最適化した。
また各振動モードはそれぞれ 8 個の基底関数で展開した。
RMC+MP 計算では、全虚時間長 β = 3.0 × 10−3 [1/K] の虚時間発展シミュレーションを行った。このとき、虚時間プロパゲータには 12 次の Multi-Product 展開を用いた。また虚時間 Path の生成数は 128 本とし、各々200,000 ステップのサンプリングを行った。
【結果と考察】 Table 1 は、基準振動解析と各 QMC 計算 (VMC、RMC+MP) による H2 Oの基本振動数と対応する実験値 [8] を示す。表中の HA は基準振動解析、VMC は線形最適化法により最適化された試行波動関数を用いた変分モンテカルロ法、VCI は Yagi らが開発した Dircet VCI 法による計算値 [1] であり、最下段は各手法による実験値からの平均絶対誤差 (MAD) である。
Table 1: H2 O の基本振動数 [cm−1 ]HA による基本振動数は実験
VMC
RMC+MP VCI[1] Expt.[8]値を大きく過大評価しており、 振動モード HA
変角
1641 1596(1)1591(1)
1576
1594.59
その MAD は 137cm−1 と大きな 対称伸縮
3838 3721(1)3662(1)3647
3657.053753(1)3760
3755.97
値になった。一方、VMC 法によ 逆対称伸縮 3942 3791(1)MAD13733(1)3(1)
11—り基本振動数は一様に改善し、MAD は 33cm−1 まで改善された。これは平均場レベルにおいて PES の非調和性と振動モード間のカップリングが考慮されたためである。一方、ポテンシャルの非調和性と振動モード間のカップリングをより精密に取り込める RMC+MP 法を適用することで、基本振動数は大きく改善し、MAD は約 3cm−1 となった。また VCI による MAD の結果との比較から、RMC+MP 法は基本振動数の実験値をより定量的に再現することが分かった。他の分子の計算結果および詳細な解析結果については、ポスター発表にて報告する。
——————————————[1] K. Yagi, et al., J. Chem. Phys. 113, 1005 (2000). [2] K. Yagi, et al., J. Chem. Phys. 127, 034111 (2007). [3]A. B. McCoy, et al., J. Chem. Phys. 123, 064317 (2005). [4] S. Baroni, et al., NIC Series, 10, 75 (2000). [5] S. A.
Chin, et al., arXiv:08009.0914v2 (2009). [6] Y. Kita, et al., In preparation. [7] U. Toulouse, et al., J. Chem. Phys.
126, 084102 (2007). [8] D. F. Smith, Jr., et al., Spectrochimica Acta, 28A, 474 (1972).