4P121
低振動分光に対応した分極水モデルの開発(京大院・理)伊藤広伸、谷村吉隆

Development of a Polarizable Water Model for Low-Frequency
Vibrational Spectroscopies(Kyoto Univ. Science)

Hironobu Ito, Yoshitaka Tanimura【序】 水は自然界において最も重要な物質の一つであり、多くの特異な性質を示す。
その特異な性質の多くは、水分子が他の水分子と水素結合することにより３次元ネット
ワークを形成するために生じると考えられている。水の巨視的な性質については長きに渡り研究され、広く知られているが、分子レベルの微視的な性質についてはまだ完全に解明されてはいない。そのため、多くの科学者らにより水の微視的な性質について今も広く研究が行われている。
分子動力学法は、ある分子力場の下にニュートン方程式を解くことにより対象となる分子の動的過程を解析する手法であり、例えば赤外振動分光スペクトルやラマンスペクト
ルのような光学スペクトルを計算・解析することに優れている。分子動力学法により計算される物理量の精度は主にポテンシャル面の精度に依存しており、精度を向上させるため
ポテンシャル面を電子状態の第一原理計算から求める事も多い。しかしながら、水分子のような液体凝縮系は分子の自由度が大きいため、どうしてもＣＰＵパワーによる制限を受けてしまい、高精度な電子状態計算を実装した系で扱える原子や分子の数は、スーパーコ
ンピュータ等を使わない限りは通常高々十数個程度である。分子動力学法は、計算精度とコストとの戦いであり、要求する物理量に応じてポテンシャルの精度を妥協していく必要があるが、近年多次元振動分光法などを始め、高精度なポテンシャルモデルを要求するシ
ミュレーションが行われることが多くなり、今後電子状態の第一原理計算を準拠とした分子動力学法が必要不可欠と思われる。
【計算、結果】 近年長谷川、谷村らにより低振動領域から高振動領域までの振動分光に対応した分極水モデル (POLI2VS) が開発された [1]。この分極水モデルは、Gaussian
03 による CCSD/aug-cc-pVQZ の高精度な電子状態計算から得られた電荷密度が利用されており、水分子の構造や周りのクラスターの構造により電荷密度が依存するモデルである。水分子内の構造による電荷密度の変化は、永久電荷密度により永久四重極子まで考慮することにより表現され、水分子が周りから感じるポテンシャルや電場に応じて生じる電荷密度は、水分子中の各原子が持つ分極率から生じる誘起単極子と誘起双極子までを考慮して表現されている。計算コストを抑えるため、電子状態計算により計算された各電荷密度は次のような関数でフィッティングされている。
∑
f (r1 , r2 , θ) =cijk ∆r1i ∆r2j ∆θk(1)
ijk∆ri = ri − req , ∆θ = θ − θeq(2)ここで、ri , θ , req , θeq はそれぞれ水分子内の OH 間の長さ、HOH の角度、平衡状態
での OH 間の長さ、平衡状態での HOH 間 の角度であり、cijk はフィッティングパラメー
タである。
この水モデルは量子力学ポテンシャルを元とした古典的シミュレーションによる制限から生じるピークのブルーシフトを除いて、約 300K における実験から得られた赤外振動分
光とラマン振動分光の高振動領域の分子内振動モードをよく再現している。一方、赤外振動分光による低振動領域の分子間振動モードは実験結果 [2] とよく一致しているが (図１)、異方性ラマン振動分光による低振動領域における分子間振動モードのスペクトルにおいて、計算により得られたシグナルと Meech らにより実験的に得られたシグナル [3]との食い違いが見られた。今回我々は、異方性ラマン振動のシグナルの食い違いが見られた原因を探求し、ポテンシャルの改善を行ったので当日報告する。

図 1: 300 K における水の赤外振動分光スペクトルの計算結果と実験結果 [1]【参考文献】[1] Hasegawa, T., Tanimura, Y., J. Phys. Chem. B 2011, 115, 5545.
[2] Bertie, J. E., Lan, Z. D., Appl. Spectrosc. 1996, 50, 1047.
[3] Heisler, A., Mazur, K., Meech,S. R., J. Phys. Chem. B 2011, 115, 1863.

4P122Reptation 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=0

j=0

ここで N は虚時間プロパゲータの分割数、dτ(= β/N) は虚時間ステップを表す。式 (3) のg は |ΨT |2 に従う分布を生成する短時間プロパゲータ、ε(x) ≡ Ψ−1
T (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)2
i=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)
MAD
137
33(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).

4P123ab initio 経路積分分子動力学法を用いたマレイン酸水素の
分子内水素結合の解析（横市大院生命ナノ１、理研 AICS2）川島雪生 1,2、立川仁典 1

ab initio path integral molecular dynamics study on theintramolecular hydrogen bond of hydrogen maleate(Yokohama City University1, RIKEN AICS2) Yukio Kawashima1,2, Masanori Tachikawa1【序論】通常の水素結合よりも結合長が短く、プロトン移
動反応の障壁が低い低障壁水素結合は、生体分子におけるプロトン移動反応において重要な役割を果たすことが示唆されており、注目を集めている。マレイン酸水素H(OOC-CH=CH-COO)− の分子内水素結合(O···H···O)はプロ
トン移動反応の障壁が低く、低障壁水素結合の研究にもよく用いられている。しかしながら、プロトン移動反応の詳
細やダイナミクスについては依然明らかにされていないことが多い。例えば、非弾性中性子散乱の実験において平面構造（Figure 1(a)や(b)）以外に面外に歪んだ安定構造(Figure1(c))が存在することを示唆する結果が得られている [1]が、理論計算では未だにその構造が確認されていない。それは、低障壁水素結合のようにプロトンが重要な役割を果たす系
を正確に記述するためには核の量子揺らぎの効果を取り込
むことが必要不可欠である[2]のに対して、それを効率よく理論計算に取り込むことが難しいからである。本研究では、核の量子揺らぎの効果を効率よく取り込むことが可能な abinitio 経路積分分子動力学法を用いてマレイン酸水素の分子シミュレーションを実行し、水素結合の構造や分子全体の運動に核の量子揺らぎの効果がどのような影響を与える
か解析を行った[3]。また、本研究では特に面外に歪んだ構
造の安定性について調査を行った。

Figure 1 マレイン酸水素の構
造: a) 静的な電子状態計算における最安定構造と b) 水素
移動反応における遷移状態の構造、ならびに、c) 面外に歪んだ構造

【計算方法】経路積分分子動力学法（以下、PIMD）では、量子粒子 N 個をそれぞれ P 個の古典粒子（ビーズ）として扱うことにより、核の量子揺らぎを表現する。本研究ではマレイン酸水素の NVT アンサンブルを発生させるために Nosé-Hoover chain を熱浴として用い、温度を 300 K、ビーズ数を 16 個、タイムステップを 0.1 fs とし、100,000 ステップの分子シミュ

レーションを実行した。電子状態計算には ωB97XD/6-31G**を用いた。分子内水素結合を形成する酸素原子には diffused 関数を追加した。本研究では核の量子揺らぎの温度効果について調べるために温度を 150 K、ビーズ数を 32 個とした分子シミュレーションも実行した。さらに、同位体効果について調べるために、D(OOC-CD=CD-COO)−についても分子シミュレー
ションを行った。また、比較のため、通常の ab initio 分子動力学シミュレーション（以下、AIMD）も実行した。
【結果・考察】まず、通常の静的な電子状態計算
でプロトン移動反応の反応障壁について計算した
ところ、0.12 kcal/mol と非常に低い障壁であることが確認できた。次に、分子シミュレーションにおけるプロトン移動反応障壁を見積もった。新たに δOH = RO3H11 − RO10H11 という反応座標を定義し、その δOH と RO3O10 の二次元自由エネルギー曲面を構築した。δOH =0 は Figure 1(b)のようにプロトンが水素結合の中央に存在し、δOH の絶対値が大きいと Figure 1(a)のようにプロトンは水素結合の両端の酸素に局在化する。AIMD と PIMD の曲面をそれぞれ Figure 2(a)と(b)に示す。AIMD ではプロト
ンは水素結合の両端の酸素原子に局在化し、プロトン移動反応の障壁が存在することが分かる。
AIMD では温度揺らぎの効果が取り込まれているが、温度揺らぎの効果だけではプロトン移動反応の障壁を克服出来ない。その一方で、PIMD では
プロトンが水素結合の中央に分布することが分かる。核の量子揺らぎの効果により、プロトン移動
反応の障壁が無くなったことが分かる。核の量子揺らぎが水素結合構造に大きな影響を与えているこ

Figure 2 マレイン酸水素の分子シとが確認できた。さらに、分子の運動の解析を行っミュレーションで得られた δOH とたところ、核の量子揺らぎの効果がマレイン酸水素

ROO の二次元自由エネルギー面: a)の面外に歪んだ構造に顕著な影響を与えていること

AIMD と b) PIMDが確認できた。解析の詳細は当日報告する。
[1] F. Fillaux, N. Leygue, J. Tomkinson, A. Cousson, W. Paulus, Chem. Phys., 244, 387-403 (1999).
[2] K. Suzuki, M. Shiga, M. Tachikawa, J. Chem. Phys., 129, 144310 (2008).
[3] Y. Kawashima, M. Tachikawa, Chem. Phys. Lett., 571, 23-27 (2013).

4P124ab initio 経路積分分子動力学法を用いた OH−水二量体への原子核の量子効果の影響（横市大院・生命ナノ 1, 理研 AICS2, 台湾中央研究院 3） ○緒方勇大 1,
川島雪生 2, 高橋開人３, 立川仁典 1

Nuclear quantum effect on OH−(H2O)2 with ab initio
path integral molecular dynamics(Yokohama City Univ. 1, RIKEN AICS2, Academia Sinica3) ○Yudai Ogata1,
Yukio Kawashima2, Kaito Takahashi3, Masanori Tachikawa1【序論】 水素結合は分子内、分子間で生じる相互作用の一つであり、あらゆる化学プロセスに影響を及ぼす。水分子は水素結合ネットワークを形成して水クラスターとなることが知られており、OH−やF−を含む水クラスターについては、今でも実験と理論の両面から盛んに研究されている[1-5]。特に、OH−の水二量体クラスター(OH−(H2O)2)は二つの構造異性体が報告されているが、その存在比は明らかになっていない。一方で OH−と水分子間のイオン性水素結合は、通常の水素結合よりも短く、強いため、このような水素結合

Conformer

を理論的に解析するためには

H5

温度効果と原子核の量子効果を考慮することが必要不可欠でSide view H1
Side view
H3
O2
O1
O3

ある[4,5]。そこで本研究では、
OH−(H2O)2 に対し、温度効果と原子核の量子効果を考慮でConformer

H2
Top view

H4
Top view

きる ab initio 経路積分分子動力学(PIMD)法を用いて、二つの構造異性体の存在比に着目

Fig. 1 Schematic illustrations of OH−(H2O)2

した解析を行った。
【理論・計算条件】 経路積分法は N 個の量子原子核を N ×P 個の古典原子核（ビーズ）として扱うことで、原子核の量子性を表現する。配置生成には温度制御法として massive Nose-Hoover
chain 法を用いた分子動力学法を採用した。PIMD 法の計算条件は温度 50 K、ビーズ数 P = 96、ステップ数 230,000 steps、ステップ幅 0.1 fs である。また、原子核の量子効果をより明確にするため、原子核を古典的に扱う従来の分子動力学(CLMD)計算も行った。計算条件は温度 50 K、ビーズ数

P = 1、ステップ数 5,800,000 steps、ステップ幅 0.1 fs である。全ての電子状態は基底関数として
TZVP を用いた RI-B3LYP で評価した。

【結果・考察】 Fig. 1 に示したよう Table 1 Dihedral angles [degrees] of 𝜙! and 𝜙! , and calculatedに OH−(H2O)2 には二つの構造 relative energy [kcal/mol] in Conformers I and II obtained by
異性体が存在する。電子状態計 electronic structure calculation.
算から得られた平衡構造の二つの二面角 H2-O1-O3-O2 (𝜙! )とH4-O2-O3-O1（𝜙! ）の値、さらに構造異性体の相対エネルギーをTable 1 に示す。Table 1 から、二つの構造異性体で𝜙! は近い値を示すが、𝜙! は大きく異なる値を示すことがわかる。すなわち、H2 原子と H4 原子が Conformer I では異なる向きを、Conformer II
では同じ向きを向く。また、全エネルギーは Conformer I の方が 0.28 kcal/mol 安定である。
これらの構造異性体の存在比を評価するため、PIMD と CLMD から得られた𝜙! と𝜙! の分布を Fig.
2 (a)と(b)にそれぞれ示す。赤い実線は PIMD による量子核シミュレーション(Qm.)を、黒い破線はCLMD による古典核シミュレーション(Cl.)を、縦に伸びる青色の一点鎖線と紫色の二点鎖線はそれぞれ構造異性体(Conf. I、II )を示す。
Fig. 2 (a) から、Cl.と Qm.ともに最大ピークは平衡構造付近に位置し、150°付近までなだらかに分布をもつ。これより、H2 原子は平衡構造から二面角が大きくなる向きにのみ大きく揺らぐことがわかる。一方で Fig. 2 (b)から、Cl.と Qm.ともにピーク位置と平衡構造の値は大きく異なり、Cl.では20°付近にピークをもち、Qm.では-60~150°に非局在化した分布をもつ。これより、H4 原子はConformer I から二面角が小さくなる向きにのみ大きく揺らぎ、原子核の量子効果を考慮することでH4 原子がより非局在化することがわかる。このように、構造異性体間の構造変化に対応する二面角において、原子核の量子効果が大きな寄与をもたらすことがわかる。
二つの構造異性体の存在比や構造変化の経路、Fig. 2 に示した分布のより詳細な議論について、

Cl.
Qm.
Conf. I
Conf. II

(a)

−120

−60

0

60

120

Dihed. [degree]φ1
H2−O1−O3−O2

istribution (arbitrary

Distribution (arbitraryunit)
istribution(arbitrary

これらの二面角の二次元分布や H5 原子の位置との関係などを用いて当日ポスターにて発表する。

Cl.
Qm.
Conf. I
Conf. II(b)180 −120

−60

0

60

120

180

Dihed. [degree]φ2
H4−O2−O3−O1

Fig. 2 One-dimensional distributions with respect to (a) 𝜙! and (b) 𝜙! (Unit in degrees)[1] W. H. Robertson, et al., Science, 299, 1367 (2003). [2] S. S. Xantheas. J. Am. Chem. Soc. 177, 10373 (1995).
[3] M. Morita, et al., Phys. Chem. Chem. Phys. 15, 114 (2013). [3] M. Morita, et al., Phys. Chem. Chem. Phys.
14, 2797 (2012). [4] Y. Kawashima, et al., J. Phys. Chem. A, 24, 5205 (2013). [5] K. Suzuki, et al., J. Chem.
Phys. 129, 144310 (2008).





溶媒効果
多成分系

（横市大生命）○兼松佑典・立川仁典
Extension of multi-component molecular4P125（Yokohama City University）○KANEMA
Abstract： Multicomponent molecular orbi
溶媒効果を考慮に入れた多成分系密度汎関数法による
analysis on solvent effect. We calculated
化学シフト・同位体シフト計算の精度検証および解析(PANO) in some solvents with this extende
those of the conventional MO method and c（横浜市立大学）	
  ◯兼松佑典, 立川仁典
the contribution of both nuclear effect andcalculation of chemical shifts.
Keywords：
MC_MO, Chemical Shift, Isoto
Accuracy check and analysis on chemical shifts and isotope
shifts calculated

by multi-component density functional theory with taking into
account of solvent effect 核の量子性を考慮した多成分系分子(Yokohama City University) ◯Yusuke Kanematsu, （PCM）
Masanori計算を新たに実装し、Tachikawa
核の量子効

本研究では、分子内に短距離水素結合をの化学シフトに着目し、H/D 同位体効果【序論】	
  NMR 化学シフトは分子の構造や化学的活性に関する情報を含有する最も PANO と標準物質であるテトラメチル
重要な分子物性の一つである。特に同位体置換によって生じる化学シフトの変化分で張 で あ る MC_B3LYP 法 に よ る 構 造 最ある NMR 同位体シフトは、原子核のおかれた化学的環境をより詳細に理解するため（!(H/D)）およびその差である同位体シに有用であり、これまでに多くの分子に対してその実験値が報告されている。量子化
基底系に 6-311++G(2d,2p), 核の基底系には IEFPCM 法を、遮蔽定数計算には CSG
学計算による NMR 同位体シフトの理論的予測もこれまでにいくつか報告されている
 PANO の H/D 同位体の化学シフト（同が、未だその手法として確立したものは存在していない。そこで我々は、NMR 化学
体シフトを表１に示す。また各最適構造
シフトおよび同位体シフトの高精度計算手法の確立を目標として、核の量子効果と溶
１より、核の量子性のみを考慮した気相計
媒効果を併せて考慮できる MC_DFT-PCM（多成分系密度汎関数法[1]–分極連続体モ
により溶媒効果を取り込んだ計算では、化現していることが分かる。表２からは、
デル[2]）を開発し、NMR 化学シフト・同位体シフトの計算と解析、および計算精度
溶液中での結合距離が気相よりも伸長す
向上のための実装を行った。
合距離と化学シフトの間に相関があるこ

表１ 気相
と重水素の

【計算詳細】	
  MC_DFT-PCM により電子と同位体置換する軽/重水素を量子的に取り扱い、化学シフト・同位体シフトの計算を行った。対象とする分子には、分子内

計算値
計算値
実験

水素結合を有し、大きな同位体シフトの実験値[3]が報告されている picolinic acid N-oxide（PANO, 図１）を採表２用した。電子の汎関数には B3LYP を用いた。電子の基底系に 6-311++G(2d,2p), 核の基底系に 1s ガウス型関数を用いた。溶媒効果計算には IEFPCM 法を、遮蔽

図１ picolinic acid N-oxide(PANO)

定数計算には CSGT 法をそれぞれ採用した。

[1] M. Tachikawa, Chem. Phys. Lett., 360, 2002

【結果】	  気相、クロロホルム溶液中、およびアセトニトリル溶液中の PANO の H/D同位体の NMR 化学シフト、および一次の同位体シフト(同位体置換部位は図 1 中に表記)を表１に示す。また各最適構造における O–H/D 結合距離を表２に示す。
表１より、気相中よりも溶液中の方がより大きな化学シフトが見られ、溶媒効果により化学シフトが増大することが分かる。次に、核の量子性を考慮しない DFT 計算で

計算値
計算値

は化学シフトを過小評価するが、核の量子性を考慮に入れた MC_DFT 計算では、それぞれの溶媒における化学シフトと同位体シフトの実験値を良く再現していることがわかる。一方、先行研究で行われた２次元の部分振動のシュレディンガー方程式を解く方法（2D-SE）による計算結果では、クロロホルム溶液中での同位体シフトがア
セトニトリル溶液中での値を上回っており、定性的にも誤った傾向を示している。
MC_DFT-PCM 法を用いて核の量子効果と溶媒効果を併せて考慮することで、より精度良く NMR 化学シフト・同位体シフトを予測できることが示唆された。
表２からは、O–H 結合距離の方がより核の量子効果が小さな O–D よりも長く、核の量子効果が結合伸長に寄与することが分かる。また、気相、クロロホルム溶液、
アセトニトリル溶液と周囲環境の誘電率が高くなるにつれて結合距離が長くなり、溶媒の分極電荷もまた結合を伸長させる寄与を有することが見て取れる。さらに表２に見られる結合長の大小と表１の化学シフトの大小は相関しており、核の量子効果と溶媒効果が O-H(D)結合距離の変化を介して同位体シフトの値に影響していることを示唆している。
	 【参考文献】[1] T. Udagawa, M. Tachikawa, J. Chem. Phys., 2006, 125, 244105.
[2] M. Cossi, G. Scalmani, N. Rega, V. Barone, J. Chem. Phys. 2002, 117, 43.
[3] J. Stare, A. Jezierska, G. Ambrozic, I. J. Kosir, J. Kidric, A. Koll, J. Mavri, D.
Hadzi, J. Am. Chem. Soc. 2004, 126, 4437.

4P126
溶液中における歯車状両親媒性分子の動的特性に関する理論的研究（横市大院・生命ナノ, 東大院総合文化, 産総研）
○増子貴子，平岡秀一, 長嶋雲兵, 立川仁典

Theoretical analysis of dynamic feature for self-assembled nanocube in aqueous solution(Yokohama City Univ., The Univ. of Tokyo, AIST)○ Takako Mashiko, Shuichi Hiraoka, Umpei Nagashima, Masanori Tachikawa【序】
親水部と疎水部をもつ両親媒性分子は, 疎水部が水中で自己組織化することにより, ミセルや脂質二重層などの多様な会合体を形成することが知られている.会合の駆動力もしくは会合機構の詳細を解明することは, 生物組織などの働きを理解する上で必要不可欠である.
ヘキサフェニルベンゼン骨格に 3 つの 3-ピリジル基を置換した歯車状両親媒性分子(図 1)は，置換基 R がメチル基の 1 では 25 %含水メタノール溶媒中にて一義的に六量体箱型集合体(ナノキューブ: 16)を形成するが, R に水素を導入した 2 は集合体を形成し

図 1 歯車状両親媒性分子

ない[1,2]. 小関らは分子軌道法および分子動力学(MD)シミュレーションにより, 気相における歯車状両親媒性分子のメチル化体と水素化体の理論計算を気相中にて実行した[3,4]. その結果, 分子軌道計算による相互作用エネルギーの観点から, π-π 相互作用よりも CH-π 相互作用が分子間に与える影響が強く, 16 は 26 に比べ安定であることが分かった[3]. また, 気相での MD シミュレー
ションにより会合体の温度効果を調べた結果，温度を上げると 16，26 共に内部体積が大きくなる
が，26 は 16 よりも低温で箱形構造を維持できなくなることが明らかとなった[4]. しかしながら,溶液中におけるナノキューブの動的安定性は未だ分かっていない. そこで本研究では, MD シミュ
レーションを用いて, 水溶媒および実験環境と同一である 25 ％含水メタノール溶媒におけるナノキューブへの自己組織化機構に関して, 溶媒分子群と歯車状両親媒性分子間の相互作用の観点から解析を行った.

【計算の詳細】結晶構造解析および NMR 測定から明らかとなったナノキュ―ブ 16 の構造特性を反映した初期
構 造 を 作成した. その周囲に SP C / E 水 溶 媒 , ま た , SP C / E 水 溶 媒 と p a r m 9 9 に よ る メ
タ ノ ー ル 溶 媒 か ら な る 25 ％ 含 水 メ タ ノ ー ル 溶 媒 環 境 を 配 置 し た . 3-ピ リ ジ ル 基

と ベ ン ゼ ン 間 の 二 面 角 パ ラ メ ー タ は HF/6-31G(d)で 得 ら れ た エ ネ ル ギ ー 曲 面 に フ
ィ ッ テ ン グ し た も の を 用 い , そ の 他 の MD 力 場 は G e n e r a l i z e d A m b e r F o r c e F i e l d( G A F F ) を 用 い た . 電 荷 は RE SP 電 荷 を 使 用 し , 周 期 境 界 条 件 の も と で A M B E R9 [ 5 ]に よ り NVT 計算を実 行 し た .

【結果】
本研究における計算対象は, 図 1 に示した歯車状両親媒性分子の単量体および六量体であり,置換基 R としてメチル基(1), 水素原子(2)を用いた. 1H NMR 測定の結果から, 単量体の分子構造は
C3 対称性を示すが, 六量体における各構成分子は C1 対称である[2]. そこで, 対称性を考慮しながら六量体の構造探索を行い, その構造・電荷を用いて新たな参照構造を生成した. 得られた六量体構造から NVT 計算を実行し, 置換基の違いや溶媒による挙動の変化を議論した.
水溶媒中での 16 および 26 における, 主成分解析の最低固有値の揺らぎ(固有ベクトル)を図 2 に示す. 16 における 3-ピリジル基の三重スタッキングとメチル基の揺らぎよりも, 26 における三重ス
タッキングと水素原子の揺らぎの方が大きいことが分かった. これは, 16 における分子間でのメ
チル基間の噛み合いが 26 における水素原子間の噛み合いよりも強いためであり, 26 では 16 に比べて自己集合しにくいと考えられる.
現在, 実験と同等な 25 ％含水メタノール溶媒中での MD シミュレーションを実行中である. 同様の解析を実行し, ナノキューブと溶媒分子間の相互作用の観点から, 自己組織化機構の要因を検

図 2 主成分解析の最低振動数の揺らぎを示す.
赤で示したものが 3-ピリジル基, シアンで示したものがメチル基(水素原子). 16 よりも 26 の方が動的揺らぎは大きい.

討したい.

【参考文献】[1]	
  S.	
  Hiraoka,	
  K.	
  Hirano,	
  M.	
  Shiro,	  M	
  Shionoya,	
  J.	
  Am.	
   Chem.	
  Soc.,	
  130,	
  14368	
  (2008).	
  [2]	  S.	
  Hiraoka,	
  K.	
  Harano,	
  T.	
  Nakamura,	
  M.	
  Shiro,	
  M.	
  Shionoya,	
  Angew.	
  Chem.,	  Int.	
  Ed.,	
  48,	
  7006	
  (2009).	
  [3]	  J.	
  Koseki,	
  Y.	
  Kita,	
  S.	
  Hiraoka,	
  U.	
  Nagashima,	
  M.	
  Tachikawa,	
  Theor.	
  Chem.	
  Acc.,	
  130,	
  1055	
  (2011).	
  [4]	  J.	
  Koseki,	
  Y.	
  Kita,	  S.	
  Hiraoka,	
  U.	
  Nagashima,	
  M.	
  Tachikawa,	
  Int.	
  J.	
  Quantum	  Chem.,	
  10,	
  1002	  (2012).	
  [5]	  D.	
  A.	
  Case,	
  T.	
   A.	
   Darden,	
   T.	
   E.	
   Cheatham,	   III,	
   et	   al.,	
   AMBER9,	
   9th	   ed.,	
   University	
   of	
   California:	
   San	   Francisco,	
  (2006).	
  4P127双直交基底第二量子化に基づく波動関数解析手法(京大院工 ∗ ，京都大学 触媒・電池元素戦略ユニット ∗∗ )○枡井 基典 ∗ ，佐藤 啓文 ∗,∗∗A Wave Function Analysis
Based on Biorthogonal Second Quantumization( Dept. Molecular Engineering, Graduate School of Engineering, Kyoto University∗ ,ESICB, Kyoto University∗∗ )○ Motonori Masui∗ , Hirofumi Sato∗,∗∗【緒言】
分子の反応性などを考える上で、化学結合が持つ共鳴構造は非常に重要な役割を果たしている。
近年、Ikeda らは分子軌道に対して結合の共鳴構造を射影する演算子を提案し、化学的直観に則した共鳴構造の重みが得られることを報告している。[1], [2] 本研究では、この方法を参考にし、対象とする軌道数や電子数の一般化を目指して研究を行った。
【理論】
本研究では、(1) 式で定義される恒等演算子を導入した。ここで χˆa σa + , φˆb σb− , · · · は双直交基底(原子軌道) 生成・消滅演算子を表し、σa , σb , · · · はスピン状態 (α , β) を表す。

1̂ =

∑

a,b,c,d Ω̂

(1)Ω̂ = 4!1 {ω̂(α, β, α, β) + ω̂(β, α, β, α) + ω̂(α, β, β, α) + ω̂(β, α, α, β) + ω̂(α, α, β, β) + ω̂(β, β, α, α)} (2)ω̂(σa , σb , σc , σd ) ≡ χ̂σd d + χ̂σc c + χ̂σb b + χ̂σa a + φˆa σa − φˆb σb − φˆc σc − φˆd σd −(3)

この恒等演算子を結合に局在化した分子軌道 ϕi , ϕ j の規格化条件の式に挟むことで、共鳴状態の射影を試みた。
ββ

1 = ⟨ϕ j ϕαj ϕi ϕαi |

∑

β

β

Ω̂|ϕαi ϕi ϕαj ϕ j ⟩(4)a,b,c,d

1 1 ∑= 5{(P i S)ab (P i S)ba (P j S)cd (P j S)dc
2 4!
4!

a,b,c,d+ (P i S)ab (P j S)ba (P i S)cd (P j S)dc − 2(P i S)ab (P i S)bc (P j S)cd (P j S)da }(5)

ここで、P i , P j はそれぞれ分子軌道 ϕi , ϕ j の密度を表す行列で、S は重なり行列を表す。

(5) 式の対称性が高いことは注目に値する。また、b, c, d について和をとると 2 つの分子軌道 ϕi , ϕ jに対する Mulliken 電荷の表式に帰着する。
1∑{(P i S)aa + (P j S)aa } = 1
4 a

(6)

【計算手法】
水分子と trans-ブタジエン分子に対して計算を行った。RHF 法から得られた正準軌道に対して

Pipek-Mezey の局在化を行い、それぞれの分子に対して O-H 結合とπ軌道に対応する 2 つの価電子軌道を得た。得られた軌道に理論を適用し、共鳴構造の重みを計算した。基底関数には TZP(5d)を用いた。
【結果・考察】Table 1 によると酸素原子と 2 個の水素原子が結合してTable 1: 水分子の共鳴構造 (一部) の重み (%)

いる共鳴構造 1 よりも酸素原子に負電荷がある共鳴構造 2

No.

の方が 2 倍ほど大きい重みを有していており、Ikeda らの

1

方法と同様に、化学的直観に則した結果を与えている。共2

鳴構造 3 や 6 のように、酸3

Present

Ikeda et al.[2]

21.1

20.6

38.5

38.5

11.6

11.2

17.8

17.8

No.

5

6

7

Present

Ikeda et al.[2]

10.8

10.5

1.6

1.5

-0.4

0.1

素が正電荷を帯びる共鳴構造もある程度の重みを有しており、このような共鳴構造を考

4

慮することにも意義があることが分かる。また、共鳴構造

Table 2: trans-ブタジエンの共鳴構造 (一部) の重み (%)

7 や Table 1 に示していない
共鳴構造の多くはわずかな負No.

Present

Ikeda et al.[2]

23.3

22.9

0.37

2.0

No.

の重みを示した。

trans-ブタジエンに関して

1

Present

Ikeda et al.[2]

22.2

21.8

24.3

24.0

3

は (Table 2)、共鳴構造 3 と 4がほぼ同じ重みを持つことか

2

4

ら、端の炭素とその隣の炭素はほぼ等しい割合で電子を保有していると言える。
総じて、本研究で用いた演算子は Ikeda らの方法と同様に化学的直観に則した結果を与え、かつ(5) 式のように高い対称性を有する式を与えた。式の対称性が高いことから、今後 3 軌道 6 電子系などへ理論を拡張することが容易であると考えられる。
【参考文献】[1] A. Ikeda, Y. Nakao, H.Sato, S. Sakaki, J. Phys. Chem. A., 110, 9028(2006).

[2] A. Ikeda, Y. Nakao, H.Sato, S. Sakaki, Chem. Phys. Lett., 505, 148(2011).
【謝辞】
本研究では Ikeda らによるプログラムを参考にしたものを使用している。

4P128
フェロセニウム系イオン液体の構造に関する理論的研究(京大院工 1、京都大学 触媒・電池元素戦略ユニット 2)○野口純樹 1, 佐藤啓文 1,2Theoretical study on the structure of ferrocenium-based ionic liquids(Dept. Molecular Engineering, Graduate School of Engineering, Kyoto University1,
ESICB, Kyoto University2)○Junki Noguchi1, Hirofumi Sato1,2【序】
イオン液体は室温で液体として存在する塩を指し、一般に耐熱性、不揮発性、導電性などの有機溶媒にはない多彩な特徴を持つことが知られている。中でも、フェロセニウム系イオン液体はFe(III)に由来する特徴的な磁性を持つことが確認されており、磁気メモリ等への応用が期待されている[1]。Fig. 1のように、弱磁場中でイオン液体を冷却して凝固させた場合は χT の上昇が小さいのに対し、強磁場中では χT が大きく上昇する。弱磁場中の凝固では多結晶が生成し、強磁場中では配向の整った単結晶が得られることが確かめられており、フェロセニウムカチオンが磁気異方性を持っているために強磁場中においては磁場の向きに沿って分子の配向が揃うと考えられている。本研究ではその異方性に着目し、理論計算から融点付近のイオン液体の液体構造の解Fig. 1.
[butyloctamethylferrocenium][TFSA]の温度磁化率曲線 (融点 299K)[1]析を行う。
【計算方法】[Fe(Cp)2]+、[TFSA]-についてそれぞれ孤立系で構造最適化を行った。基底関数は DZVP、計算方法は B3LYP 法を用い、Gaussian 09 で計算を行った。[Fe(Cp)2]+に関しては GAMESS を用いて
CASSCF 計算を行い、Spin-Orbit Coupling と 1.0T の外部磁場の影響を加えた波動関数を得た。
また、RISM-VV 計算を行い、液体中での各イオンの動径分布を求めた。それぞれの CH と CF3
は united atom として扱い、温度は 430K(融点 405K[2])、数密度は 2.183×10-3Å-3 とし、各サイトの
電荷は MK 電荷を用いた。
【結果と考察】[Fe(Cp)2]+、[TFSA]-の最適構造
を Fig. 2 に示す。[Fe(Cp)2]+は 2
つの Cp 環の配座について重なり

1.69Å

z軸

型とねじれ型の 2 通りが考えられるが、重なり型が最安定構造となり、ねじれ型は 0.4 kcal/mol 不安定であった。

[Fe(Cp)2]+[TFSA]Fig. 2.各イオンの構造また、[Fe(Cp)2]+ の最高占有軌道では Fe の dxy と dx2-y2 の 2 軌道が縮退しており、3 電子が占有する 2Eg 状態をとっていることが分かった。
次に、[Fe(Cp)2]+の CASSCF 計算を行った。Active space は 2 軌道 3 電子とし、擬縮退した 2 状
態を状態平均して計算した。外部磁場を i 軸方向(i = x, y, z)から印加したときの最安定状態の相対エネルギーΔEi と応答するスピンの期待値〈Si〉を Table 1 に示す。z 軸方向から磁場がかかった時が最も安定化が大きく、この結果から z 軸が磁化容易軸でTable 1.
あることが分かった。これは磁化率測定の実験で得られた 外部磁場に対するカチオンの応答
磁場の方向 iΔEi (cm-1)〈Si〉
容易軸方向と一致している。また、x, y 軸方向の磁場にス
x
0.00
0.00
ピンが全く応答しないことからスピンは外部磁場の方向
y
0.00
0.00
に因らず容易軸方向に向いていることが明らかになった。
z
-0.46
0.50
Table 2 には RISM-VV 計算に用いた各サイトの電荷を示し、計算の結果から得られた Fe-TFSA 各サイト間とCH-TFSA 各サイト間の動径分布関数をそれぞれ Fig. 3 とFig. 4 に示した。グラフから Fe、CH 共に最も近接しているのは O 原子であることが分かった。両者を照らし合わせると、Fe から 3.5Åのショルダーの部分では x, y 軸方向
から、Fe から 5Åのピークの部分では z 軸方向から O 原子

Table 2. 各サイトの電荷 δ
サイト
δ/|e|
Fe
0.29
CH
0.07
N
-0.34[3]S
0.78[3]O-0.50[3]CF3-0.12[3]が接近していると推測することができる。
Fe-N1

Fe-S2

CH-N1

CH-S2

Fe-O4

Fe-O5

CH-O4

CH-O5

CH-CF 3Fe-CF 3
2.5

2.0
1.5

ｇ（ｒ）

ｇ（ｒ）

2.0
1.5
1.0

1.0
0.5

0.5
0.0

0.0
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

r [Å]

r [Å]

Fig. 3. Fe から TFSA 各サイトの動径分布

Fig. 4. CH から TFSA 各サイトの動径分布【文献】[1] Y. Funasako, T. Mochida, T. Inagaki, T. Sakurai, H. Ohta, K. Furukawa, T. Nakamura, Chem. Commun.,
2011, 47, 4475 – 4477.
[2] T. Inagaki, T. Mochida, M. Takahashi, C. Kanadani, T. Saito, D. Kuwahara, Chem. Eur. J. 2012, 18,6795 – 6804.
[3] X. Zhong, Z. Liu, D. Cao, J. Phys. Chem. B 2011, 115, 10027–10040.

4P129三次元 SSSV 方程式の開発(京大院工 ∗ ，京都大学 触媒・電池元素戦略ユニット ∗∗ )○笠原 健人 ∗ ，佐藤 啓文 ∗,∗∗Development of three dimensional SSSV equation(Dept. Molecular Engineering, Graduate School of Engineering, Kyoto University∗ , ESICB,
Kyoto University∗∗ )○ Kento Kasahara∗ , Hirofumi Sato∗,∗∗【緒言】
溶液中における分子の拡散過程は，非常に重要な過程である。例えば，リチウムイオン二次電池の電解質溶液中におけるリチウムイオンの拡散は電池の性能に密接に関連している。拡散過程に対する理論的アプローチの一つとして，分子性液体の拡散方程式理論である SSSV (Site-Site

Smoluchowski-Vlasov) 方程式があり [1]，液体構造の時空間変化を追跡することが可能である。SSSV
方程式理論は短時間領域でのダイナミクスの記述には適してはいないものの，時空間変化を解析的に取り扱っているために，長時間領域におけるダイナミクスの記述において有効である。別のアプローチとして分子動力学法があり，この方法では短時間領域でのダイナミクスの良い描像を与えるが，長時間領域に関しては計算負荷のために適用が困難である。従って，SSSV 方程式理論
は分子動力学法と相補的な関係にあり，長時間領域のダイナミクスには SSSV 方程式理論が適していると考えられる。しかし，この理論は一次元溶媒和理論に基づいた定式化が行われており，一次元に縮約した時空間変化しか取り扱うことが出来ないために，電極界面近傍や外部電場下などの異方性のある系に適用することは困難である。従って，これらの系に適用するためには，まず三次元の液体構造の時空間変化を取り扱う必要がある。そこで，本研究では SSSV 方程式を三次
元溶媒和理論と組み合わせることで，SSSV 方程式の三次元化を行った。
【理論】

SSSV 方程式は溶液中での密度ゆらぎ，δρα (r) =∑ ()δ r − rαa (t) − ρα(1)aの時間発展を記述する方程式である。rαa は a 番目の分子の α サイトの時刻 t = 0 での位置，ρα は

α サイトの数密度である。SSSV 方程式を用いると，Fourier 空間 (k) における液体構造の時空間変化を記述する van-Hove 相関関数に関する微分方程式が得られる。

∂
G (k, t) = −k2 D · ρ · Φ (k) · G (k, t)∂t(2)

G(k, t) は van-Hove 相関関数 Gαβ (r, t) = ρ−1β ⟨δρα (r, t)δρβ (r = 0, t = 0)⟩ の Fourier 変換を要素とする
行列，D，ρ は各サイトの拡散係数および数密度を要素とする対角行列であり，Φ(k) は密度-密度相関関数 (δρ, δρ) の逆行列を表し，これは RISM 理論により求めることが可能である。Laplace 変換を施すことにより，式 (2) は次式のように表すことが出来る。

{}−1
G (k, s) = s1 + k2 D · ρ · Φ (k)· G (k, t = 0)

(3)

Iida らは，式 (3) 中の行列 A(k) ≡ s1 + k2 D · ρ · Φ (k) を溶質部分 (u) と溶媒部分 (v) にブロック行列として分割することで，SSSV 方程式を純溶媒系から溶質-溶媒系に拡張した [2]。本研究においても，同様の手法を適用することが出来る。まず，A(k) をブロック対角行列に分割し，さらに溶質
分子の濃度を ρu → 0 とすると，この行列は次式で表される。
()
s + k 2 Du
0A(k) =−k2 Dv · ρv · cvu (k) s1 + k2 Dv · ρv · Φvv (k)(4)

cvu (k) は溶質-溶媒間の直接相関関数である。溶媒-溶媒部分のブロック行列を，
k2 Dv · ρv · Φvv (k) = Pv (k) · λv (k) · Pv −1 (k)のように対角化し，これを用いることで A(k) の逆行列を次式で表すことが出来る。


(
)−1
s + k 2 Du0 v
 P (k) · (s1 + λv (k))−1 · Pv −1 (k)−1A(k) = 

−1
v
v
v
−1
·k2 D · ρv · cvu (k)P (k) · (s1 + λv (k)) · P (k) 


(
)−1
2
u
· s+k D(5)

(6)

式 (6) を式 (3) に代入することで，Laplace 逆変換を行うことで，Fourier 空間上における溶質のself-part に関する van Hove 相関関数 Guu (k, t) および溶質-溶媒間の van Hove 相関関数 Gvu (k, t) の解析的な表式を得る。

[]
Guu (k, t) = exp −k2 Du t[]
Gvu (k, t) = Pv (k) · Bvu (k) + Pv (k) · exp −λv (k) t · Pv−1 (k) · ρv · hvu (k)ここで，Bvu (k) は，

(
[
]
[
]) { v −1
(
)}
1


exp −λvα (k) t − exp −k2 Du t
P (k) · k2 Dv · ρv · cvu (k) α

2
u
v


 k D − λα (k)vuBα (k) =(λvα (k) , k2 Du )




[
]{
(
)}

t exp −λv (k) t Pv −1 (k) · k2 Dv · ρv · cvu (k)(λv (k) = k2 Du )
α
α α

(7)

(8)

を要素とする行列であり，hvu (k) は全相関関数である。式 (7) を Fourier 逆変換することにより，実空間上における van Hove 相関関数が得られる。式 (7) 中に含まれる関数 hvu (k) および cvu (k) は，

3D-RISM 法 [3] や MC-MOZ 法 [4] などの三次元溶媒和理論を用いることで評価出来るため，この式から三次元の溶媒和構造の時空間変化を追跡することが可能である。
実際の計算結果は，当日報告する予定である。
【参考文献】[1] F. Hirata, J. Chem. Phys., 96, 4619 (1992).
[2] K. Iida and H. Sato, J. Chem. Phys., 137, 034506 (2012).
[3] A. Kovalenko and F. Hirata, Chem. Phys. Lett., 290, 237 (1998).
[4] D. Yokogawa, H. Sato, T. Imai and S. Sakaki, J. Chem. Phys., 130, 064111 (2009).

4P130
溶液中における直鎖状分子の構造揺らぎのための
積分方程式理論（京大院工*，京都大学 触媒・電池元素戦略ユニット**）○松村 祥宏*，佐藤 啓文*,**An integral equation theory for the structural fluctuation of chain moleculesin the solution phase（Kyoto Univ.*, ESICB, Kyoto Univ.**）○Yoshihiro Matsumura*, Hirofumi Sato*,**【序論】
蛋白質のような柔軟な分子の構造揺らぎは溶媒環境に大きく影響を受ける。柔軟な分子は比較的自由に回転しうる炭素鎖を有し, そのコンフォメーション平衡（構造揺らぎ）に対す
る溶媒和効果は溶液中における化学・生物学現象を理解するための基礎的な知見である。これを調べるために, 直鎖状分子を用いた研究が行われている[1]。
本研究では溶液中の分子の構造揺らぎを調べるため, 積分方程式理論[2,3]とモンテカルロ
法を組み合わせた新規理論の開発を行った。溶液中の分子シミュレーションを行う際にしばしば生じるサンプリングの問題を改善した方法である。また, 分子シミュレーション法と同様のポテンシャルを用いるため, 水素結合や疎水性相互作用といった分子レベルの局所的な相互作用を取り入れることができる。
本方法を用いて直鎖状ブタンの構造揺らぎに対する水和効果を検討した。文献値[1]との比較検討を行い, 平衡定数や溶媒和自由エネルギーの温度依存性に関して良好な結果を得た。
本研究では, さらに, 水中における直鎖状ブタン(n-C4)から直鎖状デカン(n-C10)に至る直鎖状アルカンに理論を応用した。直鎖状分子の鎖を長くしていき, コンフォメーション平衡とそれ対する水和効果を検討した。
【理論】
溶質分子の構造揺らぎは分子内相関関数 ω で記述する。これ
は溶質の原子間距離 r に関する確率分布関数を表している。

dr dr  (| r  r | r )e
 (r;  )   dr dr e
U

V

a

ab

U

 V

b

V  V

ここで, 下付き文字は溶質分子の各原子サイトを表し（図 1）, rU ,rV はそれぞれ溶質と溶媒の座標を表す。β=1/kBT であり, kB はボ
図 1. 原子ラベルの例
ルツマン定数である。また, 系のポテンシャルは Vλ で表す。λ は摂動パラメータであり, λ=1 のときは V0=0 であり, λ=1 で完全なポテンシャル V1 となる。
溶媒和効果を考慮した分子内相関関数 ω1 を孤立分子の項 ω0 と溶媒和による変調の項 Δω に分割する, ω1=ω0+Δω。ω0 は孤立分子のモンテカルロ計算から求め, Δω は以下の OZ-type の方
程式（第一式）と closure 関係式（第二式）からなる新規積分方程式で求める。closure としては partial linear HNC-like の方程式を用いている。

ab (r;0)~ab (r )    dr'  cd (r ' ;1)  cd ,ab (r ' , r;0)cd

~ (r )  ln(1  
~ (r ))   solv (r ) for ~ (r )  0

ab
ab
ab
 ab (r;1)   ab solv~
for ab (r )  0
 ab (r )~  Δω ω であり, γ は分子内相関関数の直接項である。χ は分子内密度・密度相関
ここで, Δω0関数であり, 孤立分子のモンテカルロ計算から求める。υsolv は溶媒和による有効分子内ポテンシャルであり, RISM/HNC 方程式から求める。RISM/HNC 方程式は ω1 に依存するため, ω1に対する積分方程式理論と組み合わせて, 自己無撞着に解く。
本理論では, 溶質の構造揺らぎを顕に考慮した溶媒和自由エネルギーΔF を解析的に導出でき, 次の表式が得られる。
1
1
1
F 
 drab (r;0){2 (~ab (r)) 2 (~ab )   ab (r;  )  2 ~ab (r) ab (r;  )}2 a ,b UV  1   V  dr{caS(r ;  ) 
a ,S

1 UV
(haS (r;  )) 2 }2

これは積分方程式の結果として得られる各関数を用いて表されており, 容易に計算できる。
また, ΔF が解析的な表式で与えられているため, それぞれが物理的な意味に対応する各項に分割できる。
F  TSconf   intra   solv   mb   others
各項は構造エントロピー項(-TSconf), 分子間・分子内相互作用項(ΔΦintra), 溶媒和ポテンシャル
項(ΔΦsolv), 多体効果の項(ΔΦmb)である。最後の項(ΔΦothers)はほとんどの場合, 無視できるほど小さい。
【計算条件】
溶質としては all-atom model の直鎖状ブタン[1]と united-atom model の直鎖状ブタン(n-C4)からデカン(n-C10)までの計 8 通りのモデル分子を用いている。溶媒の水は TIP3P モデルを用
い, 温度 298.15K, 密度 0.9970gcm-3 の熱力学条件で計算を行っている。
0.03
温度 298.15K

水中（密度, 0.9970gcm-3）
孤立系

ω(r)

0.02

C4

C1

C4

C1

0.01

3.84Å

2.92Å

G

T

0.00
2.5

3.0

3.5

4.0

r/Å
図 2. 直鎖状ブタン（all-atom model）のC1-C4 間分子内相関関数
0.005

点線 孤立系
実線 水中

0.004

ω(r)

【結果と考察】図 2 には直鎖状ブタン（all-atom model）の C1-C4間分子内相関関数の計算結果を示している。2.9Å
付近にはゴーシュ(G)に対応するピークが, 3.8Å 付
近にはトランス(T)に対応するピークが存在する。
ゴーシュのピークの方がよりブロードであり, 構
造揺らぎが大きいことを示している。水中では気相中と比べてゴーシュのピークが高くなり, その割合が増加する。これは排除体積効果などの疎水性相互作用によるものと考えられる。
図 3 には n-C4 から n-C6（united-atom model）の両端間距離の分子内相関関数の計算結果を示す。
孤立系の結果と水中の結果を比較する。n-C4 から
n-C6 の分子はいずれも縮まった構造が水和により増加する。また, ここでは示していないが, n-C7 から n-C10 の場合も同様の傾向が確かめられる。分子内相関関数の計算結果から分子の平均両端間距離を計算すると, n-C4 から n-C10 のいずれの分子においても水和効果で 0.20~0.30 Å 程度縮むという結果が得られている。
溶媒和自由エネルギーの計算値 ΔF の温度・密度依存性と各熱力学状態におけるエネルギー分割の
詳細は当日に議論する。

0.003

n-C4

n-C5

n-C6

0.002
0.001
0.000
2.5

3.0

3.5

4.0

4.5

5.0

5.5

6.0

6.5

r/Å
図 3. n-C4 から n-C6（united-atom model）の両端間距離の分子内相関関数【参考文献】[1] Q. Cui, and V. H. Smith, Jr. J. Phys. Chem. B, 106, 6554(2002).
[2] T. Munakata, S. Yoshida, and F. Hirata, Phys. Rev. E, 123, 3687(1996).
[3] D. Yokogawa, H. Sato, and S. Sakaki, Chem. Phys. Lett. 487, 241(2010).

4P131

溶液中におけるピリジンの超高速緩和過程についての理論的研究(京大院理 1 , ウィスコンシン大 2 ) 中野勝博1 , 森俊文 2 , 林重彦 1 , 加藤重樹 1
Theoretical study on ultrafast decay processes of pyridine
in solution(Kyoto Univ.1 , Univ. Wisconsin2 )Katsuhiro Nakano1 , Toshifumi Mori2 , Shigehiko Hayashi1 , Shigeki Kato1
現在我々が注目している化学過程の一つとして、ピコ〜サブピコ秒のオーダーで起こる分子の電子励起状態からの無輻射緩和過程がある。これは例えば DNA の光安定性の要因とみなされるなど生物学的にも重要な現象である。また通常の化学過程は溶液中でみられるものが多く、そこでは溶媒効果により気相中とは異なる結果が生じる事が多い。本研究の目標は
このような光励起分子の無輻射緩和過程のメカニズムやその溶媒効果を明らかにすることである。
この超高速の無輻射緩和過程を理論的に説明する上で重要となるのが円錐交差 (CI) である。CI とは異なる断熱電子状態間のエネルギー差がゼロになる部分で、非断熱カップリン
グの値も非常に大きくそこを経由することで超高速緩和を実現している。CI は通常シームを形成しており、特にこのような緩和過程を特徴づけているといわれるのがエネルギーが最小となる CI(MECI) である。また溶液中では溶媒の自由度が大きすぎるため、ポテン
シャル面の代わりにそれを縮約した自由エネルギー面で考える必要がある。ここでは溶質
分子の各サイト上にかかる静電ポテンシャルを溶媒和座標とした線形応答自由エネルギー(LRFE) を用いた。我々はこの LRFE 面上で MECI を探索する方法を以前開発し、そこでは電子状態計算には MS-CASPT2 法を、溶媒和効果の計算には RISM-SCF を用いた。[1]MS-CASPT2 法は静的と動的両方の電子相関効果を含む方法で、これらの効果を考慮する
ことは二つ以上の電子状態について高精度で計算するのに必須である。また RISM-SCF により QM/MM 法ほどコストをかけずに溶質-溶媒間の相互作用を記述できる。

本研究が対象とする系はピリジン (図 1) の 1 ππ ∗ 状態からの無輻射緩和である。ピリジンは最も単純なヘテロ芳香族分子でありながら、振電相互作用やスピン軌道相互作用の効果が大きく様々な緩和経路が考えられる。またベンゼンと違い 1 ππ ∗ 状態の近くに 1 nπ ∗ 状態があることも緩和過程の複雑さの一因となっている。
過去の実験ではまず気相中のスペクトルが観測され S1 (1 nπ ∗ ) 状
態の垂直励起エネルギーは 4.74 eV、S2 (1 ππ ∗ ) 状態の垂直励起エネルギーは 4.99eV であった [2]。また系間交差収率についても

図 1 ピリジン

S1 状態に励起した場合は 0.5、S2 状態に励起した場合は 0.02 と実験的に確認された。さらに溶液中において 1 ππ ∗ 状態へ遷移した後、二成分型の緩和が観測されこの時次の五つの事実が確認された。

• 緩和の第一成分はアセトニトリル溶液中と水溶液中ではほぼ同じ速さである。
• 緩和の第二成分はアセトニトリル溶液中より水溶液中の方が速い。
• 緩和の第一成分の割合はアセトニトリル溶液中の方が大きい。
• ポンプ光の波長を 1 nπ ∗ 遷移に相当するものに変化させると、緩和の第一成分の割合が減少する。

• 溶質分子を 1 nπ ∗ 状態のないピリジニウムイオンにすると緩和の第二成分が消失する。
過去の理論ではこれらの実験事実も踏まえ次のような緩和経路が提案された [3]。

• 緩和の第一成分ははじめ S2 状態の平衡構造へと緩和しそこから遷移状態、MECI を経由してプレフルベン型へと異性化する経路。

• 緩和の第二成分はまず MECI を経由して S1 状態に遷移した後、T1 (3 ππ ∗ ) 状態への系間交差、または S0 状態への内部転換や上記のプレフルベン型以外へと異性化する経路。
このうち前者については CASSCF を用いて確認されたが、後者については具体的な計算はされなかった。そこで我々はまず参照として気相中での緩和過程を調べた。その結果、前者の速い緩和の経路については以前の研究結果とほぼ同じ結論となった。ただしプレフルベン
型より安定なヒュッケル型構造がありそちらへ異性化する経路が計算結果から示唆された。
また後者の遅い緩和の経路については S1 状態から T1 状態への系間交差が主な経路で、S0状態への内部転換はほぼ起こり得ないことが分かった。
今回は溶液中でのピリジンの、特に緩和の第二成分に相当する緩和経路での溶媒効果について調べる。手法としては上記の LRFE 面上で平衡構造や MECI などの構造最適化を行い、その自由エネルギープロファイルから緩和機構を説明する。当日は溶液中の計算結果をもとに議論を行う予定である。

参考文献[1] T. Mori, K. Nakano and S. Kato, J. Chem. Phys., 133, 064107, (2010)[2] Z. L. Cai and J. R. Reimers, J. Phys. Chem. A 104, 8389-8408, (2000)[3] M. Chachisvilis and A. H. Zewail, J. Phys. Chem. A, 103, 7408-7418, (1999)

4P132

ESP 多重極子演算子を用いた QM/MM 計算法の開発(京大院理*、ワシントン大**、京大院工***)○井上 雄介*、小杉 貴洋**、中農
山本 武志*、林 重彦*

浩史***、Development of Multipole Electrostatic Potential Operator
for QM/MM Method(Grad. Sch. of Sci., Kyoto Univ.*, Univ. of Washington**,
Grad. Sch. of Eng., Kyoto Univ.***)
Yusuke Inoue*, Takahiro Kosugi**, Hiroshi Nakano***,Takeshi Yamamoto*, Shigehiko Hayashi*

溶液や生体高分子を含む系に対する有効なシミュレーション方法として QM/MM 法が知られている。この方法は、反応中心を量子的(QM)に扱い、それ以外の部分を分子力場(MM)を用いて古典的に扱うものであり、効率的に量子効果を取り込むことができる方法として広く用いられている。従来の QM/MM 法は、一つの MM 構造の中で QM 領域の計算を行うポテンシャルエネル
ギー曲面上のシミュレーションである。しかし、実際の化学反応などは熱揺らぎの中で進行するため、自由エネルギー曲面上で議論することが望ましい。熱揺らぎの効果は、 QM/MM ハミルトニア
ンを用いて分子動力学(MD)計算を行えば取り込むことが出来る。緩和の遅いタンパク質などを含む系では長時間の MD 計算が必要となるが、QM 計算のコストが大きいために実行は困難である。この問題を克服するものとして、平均場 QM/MM 自由エネルギー法がある。この方法は 、MM 領域の揺らぎの中で QM 計算を行うもので、QM 領域の揺らぎの情報は失われるが、MD
計算のコストは大幅に軽減される。
QM/MM 法におけるハミルトニアンは、次のように三つの項に分解することが出来る：
(1)
ここで と
はそれぞれ QM 領域、MM 領域の内部のハミルトニアンであり、はQM 領域と MM 領域の間の相互作用を表すものである。この中で、はさらに静電相互作用と、van der Waals 相互作用とに分かれる。 Amber や CHARMM などの分子力場では、
MM 領域の原子電荷を、点電荷によって表現する。この表現の下で、ハミルトニアンは次の式のよ
うになる:
(2)
平均場 QM/MM 自由エネルギー法では、式(1)のハミルトニアンを用いて MM 領域の MD 計算を実行し、その平均場の中で QM 領域の電子状態計算を行う。MD 計算中に QM 領域の状態は固定されるため、通常の QM/MM­MD 計算に比べれば計算コストは小さくなるが、式(2)の第一項の期待値である QM­MM 間静電相互作用エネルギーを得るために膨大な数のクーロン積分を求める必要があり、依然として実行は困難である。この問題は、QM 領域の原子電荷も MM 領域と同様に点電荷で近似することで回避することが出来る。この方法では、 QM 領域の波動関数
によって決まる電子密度を QM 原子核上の点電荷で置き換え、QM­MM 間の静電相互作用を、
点電荷同士の静電相互作用で近似する。これによって必要なクーロン積分の数は現実的な数にまで減らすことが出来、計算コストは大幅に小さくなる。例として、 MM 領域が 1000 原子の系で 1
ns(1 ステップ 1 fs で 100 万ステップ)の MD シミュレーションを行う場合を考えると、電子密度を

近似なしで用いた場合には 100 万×1000=10 億個のクーロン積分を求めなければならない。この数は、シミュレーション時間に比例して大きくなる。一方、最小二乗法で点電荷の値を決定する
Restrained ElectroStatic Potential(RESP)による点電荷近似を用いれば、クーロン積分は点電荷の値を求める際に必要なだけであり、数千～数万個程度 (グリッドの数)で済む。この数はシミュレーション時間が大きくなっても増えることはない。
RESP による点電荷近似を用いた平均場 QM­MM 自由エネルギー計算は、既に生体分子などに応用されており、成功を収めている。溶媒を連続体で近似する RISM­SCF 法でも、溶質­溶媒間の相互作用の記述には RESP 電荷を用いており、この方法は静電相互作用エネルギーを簡便に求められるものとして有力である。多くの場合、RESP 電荷によって十分な精度で静電相互作用を記述することが出来るが、一般には複雑な構造を持つ電子の分布を、単純な点電荷の重ね合わせで表現することは、必ずしも妥当であるとは限らない。実際、RESP 電荷を用いた計算では、
電子雲による遮蔽効果を無視しているため、近距離の静電相互作用を過大評価してしまうことが知られている。また、明らかに点電荷のつくる電場は等方的であるため、不対電子などがつくる異方的な電場を再現することは不可能である。従って、例えば水素結合の方向などは、点電荷近似を用いたシミュレーションでは正しく予測することは出来ない。
本研究では、RESP 電荷におけるこれらの問題点を、ダンピング関数と双極子・四重極子演算子を導入することによって克服した。このとき、QM­MM 間静電相互作用に対応するハミルトニア
ンは次のようになる。

ここで 、 、 はそれぞれ電荷、双極子、四重極子演算子で、は電子雲による遮蔽効果を表すダンピング関数である。これらのパラメータは、RESP の方法に倣って最小二乗法によって決定する。このハミルトニアンを用いれば、RESP 点電荷近似によるものと同程度の計算コストで、高精度の QM­MM 間静電相互作用エネルギーを得ることが出来る。さらに、効率的な平均場
QM/MM 自由エネルギー計算法である ReWeigting Free Energy (RWFE) SCF 法に多重極子演算子を導入し、いくつかの系にこの方法を適用した。