3E01
多配置 Ehrenfest 法による非断熱化学動力学（リーズ大化）○斉田 謙一郎，Dmitrii V. Shalashilin

Nonadiabatic dynamics with multiconfigurational Ehrenfest
method（University of Leeds）Kenichiro Saita，Dmitrii V. Shalashilin【 序 】この十数年の計算機のめざましい発達により、化学反応動力学シミュレーションは一般的なものとなった。光化学反応のように複数のポテンシャルエネルギー曲面が関与する非断熱化学
動力学においては、電子状態理論計算における電子励起状態の取扱いの難しさとともに、非断熱領域で生じる状態間遷移を考慮する必要がある。非断熱効果を実装した分子動力学の手法として、
半古典的 Ehrenfest（SCF）法と Surface Hopping 法が広く知られている。半古典 Ehrenfest（SCF）
法では、分子系が電子波動関数によるハミルトニアン期待値から求められる仮想的な ”平均” ポ
テンシャル曲面上を運動すると考えるが、古典的な描像で直感的に理解しがたいこととともに、断熱状態がエネルギー的に解離していく際の minor path を記述することができず、反応の分岐比を議論することはできない。一方、Surface Hopping 法では分子系は単一のポテンシャル曲面の上を運動し、非断熱領域において確率論的に状態間の遷移を取り扱うため直感的な理解がしやす
いが、遷移確率の計算に Tully の Fewest Switches のようなテクニックを要したり、波束が分岐するために非常に数多くのトラジェクトリーを走らせることが必要になる。大自由度（10 原子~）系においてポテンシャルエネルギー曲面を高精度に求めて多数のトラジェクトリーを走らせることは今日の計算機をもってしても高コストのため、従来の Ehrenfest 法を改良し、複数の凍結ガ
ウ ス 波 束 を も っ て 電 子 ・ 原 子 核 波 動 関 数 を 表 現 す る 方 法 が 、 多 配 置 Ehrenfest（Multiconfigurational Ehrenfest, MCE）法［1］である。さらに現在、on-the-fly で ab initio
電子状態計算を行いながら同時に動力学計算を行う AI-MCE 法［2］の開発を行っている。
【 計 算 手 法 】AI-MCE 法［2］において、電子波動関数に市販の電子構造計算パッケージ MOLPROにより求めた断熱ポテンシャルを利用する。本研究では状態平均多配置 SCF（SA-MCSCF）法、基底関数に(aug-)cc-pVDZ 基底系を用いた。各原子核波束は Ehrenfest の原理の下、古典力学に従うが、振幅を変分的に求めることで波動関数の量子的振舞いを記述する。なお、初期状態のサンプリングは電子基底状態の調和振動解析にもとづく Wigner 分布を用いた。
【 結 果 】 エチレン分子（H2C=CH2）の電子励起状態には複数の円錐交差点が存在するため、シンプルな分子ながら非断熱動力学計算のよいベンチマークである。現在 MOLPRO 2012.1 へ実装され注目されている ab initio multiple spawning（AIMS）法［3］では、最初 1000 個だったトラジェクトリーが次々と分岐（spawning）するため最終的に 5000 トラジェクトリーの計算を要しているが、AI-MCE 法ではわずか 120 個のトラジェクトリーで定性的な結果が得られ、500 ト
ラジェクトリーでほぼ収束を見せた（図１）。
(a)τ = 147 fs(b)τ = 112 fs(c)

図１. エチレン分子ππ*励起後の状態の割合の時間変化．黒が S1 状態、赤が S0 状態をそれぞれ示している．(a) AIMS 法、5000 トラジェクトリー、(b)AI-MCE 法、120 トラジェクトリー、(c)AI-MCE 法、
500 および 1000 トラジェクトリー．

	  次 に 、 ヘ テ ロ 環 化 合 物 で あ る ピ ロ ー ル 分 子（C4H4NH）の光励起水素脱離反応に AI-MCE 法を適用した［4］。S1（1A2）状態の FC 領域に置かれた波束はただちに（< 50 fs）最初の円錐交差点に達し水素原
子とピロリルラジカルに解離した。また、ピロリルラ
ジカルにも複数の円錐交差点があり、水素脱離後も複雑な状態の割合変化が見られた（図２）。これらの非断熱遷移領域では面外・変角振動モードが顕著に現れており、光フラグメント分光観測結果［5］を説明できる。
これまでの動力学理論研究［6］では詳細に議論されて

図２. ピロール分子 S0→S1 励起後の状態の重みの時間変化．
赤が 1A2、
青が 1A1、緑が 1B1、また、実線が S1、破線が S0、点線が S2 状態をそれぞれ示している。

いなかった部分である。
【参考文献】[1] D. V. Shalashilin, J. Chem. Phys., 132, 244111 (2010).
[2] K. Saita, D. V. Shalashilin, J. Chem. Phys., 137, 22A506 (2012).
[3] B. G. Levine, J. D. Coe, A. M. Virshup, T. J. Martinez, Chem. Phys., 347, 3 (2008).
[4] K. Saita, M. G. D. Nix, D. V. Shalashilin, Phys. Chem. Chem. Phys., in press(DOI:10.1039/C3CP51199E).
[5] B. Cronin, M. G. D. Nix, R. H. Qadiri, M. N. R. Ashfold, Phys. Chem. Chem. Phys., 6, 5031 (2004).
[6] V. Vallet, Z. G. Lan, S. Mahapatra, A. L. Sobolewski, W. Domcke, J. Chem. Phys., 123, 144307(2005).

3E02
変分経路積分分子動力学法と電子状態計算の結合（金沢大理工*，分子研**）◯三浦 伸一*，河津	
  励*, **

Variational path integral molecular dynamics method
combined with electronic structure calculations(Kanazawa Univ.*, Institute for Molecular Science**)◯Shinichi Miura*, Tsutomu Kawatsu*, **	  量子モンテカルロ法は、多自由度系の基底状態および低励起状態を高精度に計算することができる一群の確率論的シミュレーション技法である。例えば、変分モンテカ
ルロ法や拡散モンテカルロ法と呼ばれる方法がその範疇に入る。液体ヘリウムのような量子液体に始まり、断熱ポテンシャル上での分子振動やトンネル現象、さらには原子・分子およびその集合体の電子状態もその射程に含まれつつある[1,2]。演者はこれまで、多粒子系の基底状態を高精度で計算することが可能な変分経路積分法[2]を取り上げ、その動力学的なシミュレーション手法を開発してきた[3,4,5]。この方法
を変分経路積分分子動力学法と呼ぶ。変分経路積分法は、拡散モンテカルロ法(DMC)と同程度の高精度計算が可能であり、またDMCがあまり得意としないハミルトニアンと交換しない物理量の計算も容易に実行することができる。また分子動力学法を用いることにより、ウォーカー別の並列化のみならず、各ウォーカーの高並列化も可能となり、超並列コンピュータに適した方法といえよう。本研究では、分子集合体の量子振動状態を念頭におき、断熱ポテンシャル面を非経験的分子軌道法により記述する手法を開発した。	
 
まず変分経路積分法について簡単にまとめておく。多粒子系の基底状態を試行関数	 で表す。この状態から以下の関係式より、厳密な基底状態	
 | 0 を取り出すことができる：	
 |⇥0 = lim e

2 Ĥ

|

T

	
 
ここでβは虚時間と呼ばれる量である。この基底状態の内積を考えよう：	 ⇥⇤⇥0 |⇥0 ⇥ =
=

T |e

Ĥ

|

T⇥

dRdR⇥
Ĥ

T |R⇥ R|eĤ|R⇥ ⇥ R⇥ |

T⇥

	
 

⇥

積分の中に表れている行列要素	
  R|e |R ⇥ は、逆温度βでの密度行列である。つまり有限温度での経路積分法を用いて、上記の内積を表現できることを意味している。離散化した経路積分を用いると	
 ⇥⇥0 |⇥0 ⇤ =となる。ここで	
 S({R

(s)

};

···

dR(0) · · · dR(M )

T (R

(0)

)e S({R

(s)

};

)

T (R

(M )

)

	
 

) は、離散化した虚時間経路積分表示での作用、Δτ = β / Μは虚時間ステップである。十分「長い時間」βをとれば、時間β/2 あたりの分布は厳ï1V ( cm )密な波動関数に従っていることになる。また、経路積分表示の分配関数を「古典的な高分子系」の分配関数と読み替えたように、波動関数の内積も高分子系との対応をつけることができる。分配関数の場合との違いは、高分子が開いているということと、
その高分子の両端の分布が試行関数の影響を受けていることである。このことを念頭において、有限温度での経路積分分子動力学法と同様に仮想質量、仮想運動量を導入することにより分子動力学法を構築することができる。この手法が変分経路積分分子
動力学法である。また各虚時間で核は断熱ポテンシャル面上にある。この断熱ポテン
シャルを非経験的分子軌道法と結合することにより電子状態計算から直接求める手法を第一原理変分経路積分分子動力学法と呼ぼう。	
 	
  本研究ではモデル系として水１分子を選んだ。試行波動関数はガウス型のもの4000を採用した。まず、断熱ポテンシャルとして経験的なポテンシャルを用いた場合の結果を示す。図は水分子のポテンシャ
3000

ルを虚時間の関数として示したものである。H2O、D2Oともに時間が経過するにつれて一定値に収束していく様子が見て取2000れる。これはシステムが基底状態に到達したことを示している。D2Oの場合は、H2O
00.00050.001
より基底状態に到達する時間が余分にかï1o (K )
かっているが、これは試行波動関数の変分パラメータがH2Oに対して最適化され
図：	
  虚時間の関数として水分子のポテンシ
たものを用いているためである。つまりャルを示す。上の線は H2O、下の線は D2O
貧弱な試行関数を採用しても十分なβをの結果である。
用いれば精密な結果が得られることにな
り、変分経路積分法の頑強性が示されている。一方で、電子状態計算と結合する場合
は、計算効率を最大化するためになるべく短いβで収束することが望ましい。つまり
試行波動関数の事前の最適化が重要であることを意味している。電子状態計算と結合した計算結果は当日示す。
参考文献	
 [1] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73,
33 (2001).
[2] D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
[3] S. Miura, Chem. Phys. Lett. 482, 165 (2009).
[4] S. Miura, Comp. Phys. Comm. 182, 274 (2011).
[5] S. Miura, in: S. Tanaka, S. Rothstein, W. Lester (Eds.), Advances in Quantum
Monte Carlo, 177 (2012).

3E03

第一原理経路積分インスタントン法の開発(分子研*、金沢大理工**) ○河津

励*,**、三浦伸一**

Development for the ab initio path integral instanton method(IMS**, Kanazawa Univ.**) ○Tsutomu Kawatsu*,**, Shinichi Miura**【序】水素原子等のトンネル現象は、アンモニアの反転運動やメチル基の回転運動等でよく知られている。このようなトンネル現象による状態間のカップリングはトンネ
ル分裂としてマイクロ波吸収などで直接観測されている。トンネル分裂を計算する手法の一つとして、経路積分インスタントン法[1]が知られている。経路積分の計算手法
として、離散化した虚時間スライスの分子構造等のポテンシャルを独立に計算し、時間スライス間は調和的に相互作用している離散化経路積分法が研究されており、そのポテンシャルを用いたサンプリングから分配関数を計算する経路積分分子動力学法
や経路積分モンテカルロ法が知られている。この離散化経路積分法を用いてインスタントン法を実装する手法が Althorpe らによって開発された[2]。我々は彼らの方法を踏襲しつつ、独自に離散化経路積分インスタントン法のプログラム開発と発展を行っている[3]。本研究では、量子化学計算とインスタントン法を結合した第一原理経路積分
インスタントン法を提案し、モデル計算を実施した。
【計算】インスタントン法によるトンネル分裂計
算では、Fig.1 に見られるような対称二重井戸問
題において最も理論的にクリアな表式を与える
ため、現状ではこの問題に限定した開発を行っている。テスト用のモデル系としては、この問題で最もよく知られているアンモニア傘反転運動を主に用いており、そのほかにいくつかの比較的小規模な分子系についても取り扱っている。対称二重井戸問題の解法は非対称なポテンシャルのケースには直接応用できないが、メチル基回転のような対称三重井戸問題に対しても、３つの二重井戸を仮定することで近似的に応用できる。
Fig.1. Double-well model。
我々のプログラムは MPI で虚時間スライスのポ
テンシャル計算の単位で並列化されており、各プロセスが独立に動く量子化学計算ソ
フトウェア GAMESS を呼び出して、計算を行う仕様になっている。そのため、本体のプログラムの並列度とは別に GAMESS の並列化も合わせて扱うことができる。

【結果と考察】大量の構造サンプリングが必要な分子動力学法やモンテカルロ法とは異なり、インスタントン法では一度の構造最適化のみを行う。そのため、比較的コストの高い量子化学計算の方法や大きな基底関数を用いることができる。
アンモニア系を用いた詳細なテストでは、基底状
態や遷移状態の構造および遷移状態のポテンシャルバリアを正しく再現できるレベルの量子化
学計算の方法を選ぶことで、実験値に近いトンネ
ル分裂を再現でき、逆にそれらを再現しない計算
方法の組み合わせでは不十分であることがわかった。これは二重井戸のトンネルの問題がバリアの高さとその幅に依存しているためであると考 Fig.2. One-dimensional mass-weighted
tunnel path among different atomic
えられる。
weight。
一方、トンネル現象はトンネルする原子の質量に強く依存するが、これは原子の座標を原子質量の平方根で重みづけした座標系を用いることで、上記のバリア幅の伸び縮みの一部として考えると都合が良い。この座標を用いて考えると、同じ断熱ポテンシャルを持つ同位体分子同士の違いが主にトンネ
ル距離の違いとして説明でき、その距離に伴いトンネル時間も伸び、トンネル分裂も減少することが確認できた。計算した一次元トンネル経路(インスタントン)とトンネ
ル距離 SM+1 の同位体効果による違いを Fig.2 に示した。水素を同位体置換した ND3はノーマルのアンモニア分子 NH3 比べて傘反転のトンネル距離が 25%ほど長く、トンネル分裂は 1/16 ほどになる。窒素を置換した 15NH3 でもわずかにトンネル距離が長くなり、5%ほどトンネル分裂が小さくなった。これらのトンネル分裂の間の比率は実験値をよく再現している。
[1] Mil’nikov G V and Nakamura H 2001 J. Chem. Phys. 115 6881.
[2] Richardson J O and Althorpe S C 2011 J. Chem. Phys. 134 054109.
[3] Kawatsu T and Miura S J. Phys. Conf. Ser. in press.

3E04
電荷移動状態を経由した励起エネルギー移動理論の開発と応用(神戸大院・システム情報*)○藤本和宏*

Theory of Excitation-Energy Transfer and its Application(Kobe Univ.*)○Kazuhiro J. Fujimoto*【序】 励起エネルギー移動（EET）は分子ペアや分子集合体において観測される現象であり、その特性は植物の光合成に代表されるように生体内での反応において広く用いられている。生体内での EET 機構の本質的理解は、高効率の太陽電池の開発に対しても役立つと考えられている。
EET 機構の説明は、主に Förster 理論と Dexter
理論を用いて行われきた。Förster 理論では分子間

ET1

HT1

HT2

ET2

の擬クーロン相互作用によって、Dexter 理論では
分子間における電子のスピン交換相互作用によって EET 機構が説明される。これらに対し、1994 年に Harcourt らは電荷移動状態を経由した EET という新たな機構を提案した[1]。Fig. 1 に示すように、
Harcourt 機構では、ドナー分子が励起状態となる
始状態とアクセプター分子が励起状態となる終状態以外に、各分子が電荷分離した中間状態を含めて EET が説明される。Harcourt らの報告以降、この機構に関する研究は行われてこなかったが、私

Fig. 1. Schematic illustration of EET via CT states.
Blue and red lines show HOMOs and LUMOs ofmonomer (I or J), respectively.

は Transition-Density-Fragment Interaction（TDFI）法と Transfer Integral（TI）法を組み合わせた手法（TDFI-TI 法）によって Harcourt 機構が定量的に記述できることを示してきた[2]。
【結果】 本発表では、TDFI-TI 法の説明およびテトラセン誘導体分子結晶に対する応用計算を示す[3]。TDFI 法は電子的カップリングの一つである擬クーロン相互作用を記述するために開発された手法である。これまでに、キサントロドプシンにおける EET[4]やエキシトン結合型の CD ス
ペクトル[5]に対して TDFI 法は適用されてきた。これらの研究の中で、TDFI 法は双極子・双極子
近似や TrESP 法よりも精度よく擬クーロン相互作用を記述できることを示してきた。また、分子間の self consistency を満たすために採用している DFI 法[6]が、計算精度の向上に寄与していることも示してきた[5]。従来の TDFI 法にスピン交換相互作用を加え、さらに TI 法と組み合わせたものが TDFI-TI 法である[2]。スピン交換相互作用の追加により Dexter 機構への適用を、TI 法との組
み合わせにより Harcourt 機構への適用を可能としている。昨年度の発表では電子的カップリングの摂動計算に対して TDFI-TI 法を使用したが、本研究では以下に示すハミルトニアン行列の構築に対して TDFI-TI 法を使用した。

 E1 VEET VET VHT VEET E2 VHT VET H VET VHT E3 VCT  VHT VET VCT E4 

(1)

ここで、対角成分である E1 と E2 は分子内励起のエネルギーを、E3 と E4 は電荷分離状態のエネルギーを表す。非対角成

(a)

(c)

(b)

(d)

分である VEET 、 VET 、 VHT 、 VCT は（擬クーロン＋スピン交換）項、電子移動項、ホール移動項、電荷移動項をそれぞれ示し、各項は TDFI-TI 法によって求めた。
本手法を 2 種類のテトラセン誘導体分子結晶へ適用した[3]。Fig. 2 に示す 2 種類のテトラセン誘導体（BUR と
BUY）は結晶中で異なる分子間配置をとり、それぞれ異なる吸収波長を有することが知られている。こうした現象はクリスタロクロミーと呼ばれている。式(1)を対角化して励起エネルギーを計算した結果、実験値を精度よく再現することに成功した。テトラセン誘導体の単量体から六量体に対して TDFI-TI 法を適用した結果を Table I に示す。BUR

Fig. 2. (a) BUR monomer, (b) BUY
monomer, (c) BUR hexamer, and (d) BUY
hexamer.

と BUY では多量体効果による赤方シフトの具合が異なる
ことが分かった。また、式(1)の非対角項のうち VEET だけを考慮して励起エネルギーを計算したところ（2-state model）、単量体における値と同様のものが得られた。この結果から、本系では一般的に電子的カップリングの主因と考えられる VEET よりも VET 、 VHT 、 VCT の寄与の方が大きいことが分かった。
本発表では、EET 機構以外にクリスタロクロミー機構の詳細についても述べる予定である。
Table I. TDFI-TI absorption energies of the BUR and BUY N-mers (eV).a
Monomer

Dimer

BUR

2.56

BUY

2.68

a

Trimer Tetramer Pentamer

HexamerExp.

2.37 (2.55b)

2.29

2.25

2.23

2.23 (2.23c)

2.18

2.65 (2.67b)

2.62

2.61

2.61

2.61 (2.62c)

2.56

SAC-CI data are employed. b 2-state model result. c Absorption energy calculated in crystal.

References:[1] R. D. Harcourt, G. D. Scholes, and K. P. Ghiggino, J. Chem. Phys. 101, 10521 (1994).
[2] K. J. Fujimoto, J. Chem. Phys. 137, 034101 (2012).
[3] K. J. Fujimoto and C. Kitamura, J. Chem. Phys. in press.
[4] K. J. Fujimoto and S. Hayashi, J. Am. Chem. Soc. 131, 14152 (2009).
[5] K. J. Fujimoto, J. Chem. Phys. 133, 124101 (2010).
[6] K. Fujimoto and W.-T. Yang, J. Chem. Phys. 129, 054102 (2008).

3E05
	
 
分子間の引力への軌道間相互作用の寄与：	
 水素結合、カチオン/π	
 相互作用、CH/π	
 相互作用の解析	 （産総研）○都築誠二，内丸忠文	
 Contributions of orbital-orbital interactions to the attraction between molecules:
Analysis of hydrogen bond, cation/pi and CH/pi interactions(National Institute of Advanced Industrial Science and Technology)○Seiji Tsuzuki, Tadafumi Uchimaru
	
 
【序】	
 
	
  分子の間には種々の原因による分子間力（相互作用）が働く。分子間力は交換反発力や電荷移動力のような分子軌道の重なり（軌道間相互作用）が原因の短距離力（分子軌道の重なる短距離でしか働
かない相互作用）とクーロン相互作用が原因であり、軌道の重なりが生じない遠距離でも働く遠距離力に分けることができる。分子が元々持っている電荷の間の相互作用である静電力、誘電分極による引力である誘起力、電子の相関運動を原因とする分散力は遠距離力に分類される。分子間相互作用へのそれぞれの分子間力の寄与を知ることは分子間相互作用の性質（相互作用の強さ、距離依存性、方向依存性）を理解する際に重要である。このため分子間相互作用への各分子間力の寄与の解析がこれまでにも試みられてきた。軌道計算の際に拘束条件を用いるエネルギー分割法を用いた場合、かなり大きな電荷移動エネルギーが計算されることから水素結合などの分子の相互作用において、軌道間の相互作用が引力に大きな寄与をしているとしばしば主張されてきた。今回、超分子法で計算された相互作用エネルギーと	
 distributed	
 multiple	
 を用いた解析などから得られた静電力、誘起力の寄与を比較することで軌道間の相互作用の分子間相互作用への寄与を検討したので、その結果を報告する。
	
 
【方法】	
 
	
  分子軌道法計算には Gaussian09 プログラムを使い、分子間距離を変えて相互作用エネルギーポテ
ンシャルを計算した。全相互作用エネルギー（EMP2）は MP2/aug-cc-pVTZ レベルで計算し、基底関数重ね合わせ誤差は counterpoise 法で補正した。MP2 法で計算したモノマーの電子状態から各原子の上に置く多極子（distributed multipole）を計算し、distributed multipole の間の相互作用として静電力の寄与（Ees）を計算した[1]。distributed multipole の作る電場と原子の分極率から誘起力の寄与（Eind）を計算した[2]。HF 法で計算される相互作用エネルギー表１ 超分子法で計算されるエネルギー（EHF）は表１に示すようにHF 法で計算されるエネルギー（EHF）

静電力、誘起力および軌道間

静電力の寄与（Ees）の相互作用（交換反発力、電誘起力（誘電分極）の寄与（Eind）

荷移動力）の寄与（Eshort）か

軌道間の相互作用（短距離相互作用）の寄与（Eshort）

らなるので、軌道間の相互作（交換反発力、電荷移動力）

用（短距離相互作用）の寄与

電子相関の寄与（Ecorr）

は Eshort = EHF – Ees – Eind と分散力の寄与し て 計 算 し た 。 MP2 法 と電子相関が静電力、誘起力に与える影響

HF 法で計算されたエネルギーの差である電子相関の寄与（Ecorr）には、分散力の寄与と電子相関の補正が静電力、誘起力に与える影響が含まれる。中性分子の相互作用では電子相関が静電力、誘起力に与える影響は小さいので Ecorr の大部分が分散力である。
【結果】	
 
	
  図１に計算された水二量体（φ=	 4.3˚,	
 θ=	
 116.7˚）の全相互作用エネルギー、静電力、誘起力、
電子相関、軌道間の相互作用の寄与を示す。水素結合の引力の大部分は静電力であり、誘起力、分散力も引力に寄与している。一方、軌道間の相互

作用（交換反発力と電荷移動力の和）は全体と5

しては常に反発になっている。	
 

H

H θ
O

	
  図２にリチウムイオンとベンゼンのカチオオンとベンゼンの間には水素結合と比べて非常に強い引力が働いている。カチオン/π相互-1-5

作用では引力の大部分は静電力と誘起力であ

-10
2.4

静電力や誘起力を若干大きく計算するので、カ力の大部分は分散力である。軌道間の相互作用

3.6
R/Å

-1

-20

-60

作用のいずれの場合も、軌道間の相互作用が引

MP2
es
ind
corr
short (rep + CT)
1

力の重要な原因になっていると主張されることがしばしばある。しかし、今回の解析ではい

5

2

3

H

全体としては斥力となっており、軌道間の相互

C

5

H
H
R

E /kcal mol

-1

H

0
MP2
es
ind
corr
short (rep +CT)

[1] A. J. Stone and M. Alderton, Mol. Phys., 56, 1047
(1985).
[2] A. J. Stone, Mol. Phys., 56, 1065 (1985).

4

R/Å

ずれの相互作用の場合も軌道間の相互作用は

【文献】

4.8

0

	
  水素結合、カチオン/π相互作用、CH/π相互

にくい。	
 

4.4

R

-40

は	
 CH/π相互作用の場合も斥力となっている。
	
 作用が引力に重要な寄与をしているとは考え

4

20

E /kcal mol

は静電力、誘起力の引力への寄与は小さく、引

3.2

Li+

様に全体では反発になっている。	
 用の解析結果を示す。CH/π相互作用の場合にH40いる。また、軌道間の相互作用は水素結合と同	  図３にはメタンとベンゼンの	
 CH/π相互作

2.8



チオン/π相互作用では	 Ecorr	 が反発になって

O

MP2
es
ind
corr
short (rep + CT)る。リチウムイオンが強い電場を持つため、誘起力は非常に大きい。HF	
 法は	
 MP2	
 法と比べてH

0

E /kcal mol

ン/π相互作用の解析結果を示す。リチウムイ

R

φ

-5

3

4

5
R/Å

6

3E06
実空間グリッドによる電子密度汎関数法の超並列化とその評価(東北大学) ◯高橋英明*、三木雄詩、森田明弘
Development and assessment of the massively-paralleldensity functional theory using real-space grids(Tohoku Univ.) ◯Hideaki Takahashi, Yuji Miki, Akihiro Morita【はじめに】

近年の電子状態計算や分子シミュレーションの急速な発展は、計算機の進歩に負うところが大きい。最新の計算機の速度向上は CPU 単体の高速化よりも、むしろ通信機器で接続された CPU の集積度の増大によっており、これは、「京」コンピュータをはじめとする多くのスパコンに共通する設計指針となっている。こうしたコンピュータは並列計算機とよばれ、分散メモリ上でタスクを並列処理することにより計算を高速化する。並列計算機は、現在、研究室から計算機センターに至るまで、高速演算サーバーの主流になっている。通常、並列計算では、
CPU 間の情報通信は、並列用のライブラリを介してユーザーによって能動的に制御される。高速化の成否は並列化効率に支配されるので、計算の手段としてどのような方法論やアルゴリズムを採用するかということは大規模計算の実行においては決定的に重要である。従って、電子状態計
算の方法論は、並列計算機のアーキテクチャに適合するように設計するのが望ましい。本研究では、Kohn-Sham の密度汎関数法 (KS-DFT)による電子状態計算の並列化を考える。
高い並列化効率を実現するための要諦は、CPU 間の通信量を可能な限り小さくすることにある。
KS-DFT 法では、電子間の交換や相関エネルギーを実空間上の局所、及び準局所なポテンシャルで記述するので、実空間を分割する並列化によって高い並列化効率が得られると期待される。本
研究では、実空間グリッドによる KS-DFT 法とその並列化[1]を概説し、LCAO 基底による計算結果と比較することによりその性能を評価する。また、超並列化 QM/MM シミュレーションによる自由エネルギー計算の現状を紹介する。
【方法論】

KS-DFT の基礎方程式は、1 電子波動関数とその固有値をそれぞれ、ϕ i , ε i と

すると、

(1)
で与えられる。式(1)の左辺の括弧内は左から順に、電子の運動エネルギー、Hartree ポテンシャ
ル、外部ポテンシャル、及び電子の交換相関ポテンシャルを表す。また、式(1)において、演算子
υ H と υ xc の電子密度 n の依存性を明示した。我々の実空間 KS-DFT では、運動エネルギー演算子

は典型的に 4 次の差分法で表示する。また、 υ H は Poisson 方程式を数値的に解くことによって構築する。Poisson 方程式において、 υ H のラプラシアンは運動エネルギーと同様、4 次の差分表示を用いる。原子・分子の計算において、外部ポテンシャル υ ext の主要項は原子核の引力ポテンシャ
psルである。我々は、これらをノルム保存の分離型擬ポテンシャル υ ext
によって記述している。このポテンシャルは、長距離的な局所項と短距離で減衰する非局所項からなり、実空間表示においても擬ポテンシャルの非局所性は空間の狭い領域に限られる。また、コア付近での擬ポテンシャルの急激な変化を再現するために、2 重グリッド法を採用している。式(1)左辺の括弧内の最後の項

υ xc の扱いは LCAO 法の場合と同じであり、GGA 汎関数を使用したとしても、 ∇n を考慮するだけでよい。このように、式(1)の左辺の演算子の多くが局所あるいは準局所的であることが分かる。
これは、実空間表示による KS の有効ハミルトニアンの表現行列が対角優勢になることを表しており、空間分割による並列化によって通信コストが低減できる。
実空間における空間分割のスキームを図１に示す。分割後のサブセル内のデータにおいて、通信に関わるのは他のサブセルに接する部分の波動関数や電子密度の値のみである。また、これらの通信はローカルに行われるので、大規模系においても通信量は線形に増大するのみである。また、各サブセル上の計算は OpenMP や Openacc によってスレッド並列することも可能である。
【計算の詳細】氷の構造を持つ N 個の水分子のクラスター (H2O)N について、それらの電子状態計算を MPI 並列による実空間 KS-DFT によって実行し、計算時間を計測した。また、それをスレッド並列による Gaussian09(G09)の計算時間と比較した。計算には、ノード当たり 16 コアを有する Intel Xeon (E5-2650) 4 ノード(総コア数 = 64 コア)からなる並列計算機を使用した。ただし、G09 では、1 ノード(16 コア)での計算時間を実測し、その速度が線形に増大するものと仮定して64 コアでの時間を推定した。実空間 DFT の実空間セルは立方体とし、x, y, 及び z 軸方向にそれぞれ、120 あるいは 140 個のグリッドを等間隔に配置した。グリッド幅は、典型的に 0.152 Å とした。また、BLYP を交換相関汎関数として用いた。この設定での実空間 KS-DFT は、LCAO 基底の aug-cc-pVDZ や aug-cc-pVTZ とほぼ同等の計算精度を実現するので、G09 では aug-cc-pVDZ
基底を用いた。リファレンスとして、京コンピュータの 512 ノードを駆動する計算も実行した。
【結果と考察】図 2 に N = 11〜61 の水クラスターについて、KS-DFT の SCF 1 ステップ当たりの平均の経過時間をプロットした。これらのサイズの系では軌道の直交化にかかる時間の割合が小さいために、実空間 DFT の計算時間は水クラスターのサイズに対してほぼ線形に増大する。
これに対して、LCAO 計算では、実空間法と比較して計算時間が顕著に増大する。今回の計測では、LCAO 計算は共有メモリー上のスレッド並列で行われたので、並列化に伴う通信コストは発生しない。従って、原子基底の 4 中心 2 電子積分に付随する処理が LCAO 計算のコスト増大の主要因と考えられる。電子密度を別途、基底関数で展開する手法を用いれば、この状況は緩和されるだろう。また、今回の実空間 DFT の京コンピュータへの実装では、実効性能は低いものの、極めて高速な計算が可能になることが分かった。現在、超並列化 QM/MM シミュレーションが稼働しており、反応自由エネルギー変化の計算が可能になっている。
[1]

H. Takahashi et al, Advances in Quantum Chemistry, 59, 283-351 (2010).

図１

実空間セルの空間分割

図2

SCF 1 ステップ当たりの平均の経過時間

3E07
水和における π 電子の揺らぎの寄与の自由エネルギー解析（東北大院・理）○鈴岡大樹、高橋英明*、木村慎太郎、森田明弘

Free energy analyses of the fluctuation of π electrons in a polar solvent(Tohoku Univ.) ○Daiki Suzuoka, Hideaki Takahashi*, Shintaro Kimura, Akihiro Morita

π 電子は芳香族分子や生体分子など様々な分子に遍在する基本的な電子構造の一【はじめに】つである。π 電子は空間的に大きな広がりを持ち、外場に容易に応答して分極する。π 電子はこのような性質により、CH/π 相互作用を始めとして、他のσ結合と明確な分子間相互作用を形成する。
例えば、アミノ酸の芳香環とσ電子との相互作用は、タンパク質のフォールディングにおいて重要な役割を担うと考えられている。また、芳香環を持つ最も単純な分子であるベンゼン分子は、代表的な無極性分子であるにも拘らず負の水和自由エネルギーを持つことが知られている。これは、定性的には極性分子である水分子の OH 結合により瞬間的に誘起されるベンゼン分子中の π 電子の揺らぎによる安定化の為であると予測されている。しかし、これまでの所、凝縮系の熱力学量におけるσ-π 相互作用の寄与についての定量的な研究は行われていない。
本研究は、ベンゼンの水和自由エネルギーにおける π 電子の揺らぎの寄与を定量的に解析することを目的とする。π 電子の揺らぎを扱うためには量子力学に基づく手法が不可欠である。そこ
で、摂動論に基づく QM/MM 法と溶液論を結合し、溶媒和自由エネルギーを π 電子とσ電子のそれぞれの寄与に容易に分割することが出来る新たな方法を開発した。
【理論と方法】

著者らは、溶媒和自由エネルギーの計算の方法論として、QM/MM 法とエネル

ギー表示の理論（ER 法）[1]を結合した QM/MM-ER 法を開発した[2, 3]。QM/MM-ER 法の枠組みの中では、溶質の電子密度が固定された中間状態を導入することにより、溶媒和自由エネルギーΔµを二体的な相互作用による寄与 Δµ と電子密度が揺らぐことによる寄与δµに分割する。すなわ
ち Δµ = Δµ + δµ 。Δµ は電子密度を固定した溶質が溶媒中に移動する際の溶媒和自由エネルギーに相当する。一方、δµは溶質の電子密度が、固定された電子密度から溶媒の運動に応答して分極することにより生じる自由エネルギー変化を示す。
本研究では、二次までの摂動論によりδµを π 電子とσ電子の寄与に分割する。気相中のハミル
トニアン H0 を無摂動ハミルトニアンとし、溶媒分子が作る電場 Vpc を摂動とすると、電子分極に起因するエネルギーは二次の摂動エネルギーで表される。無摂動系の一電子波動関数 {φi } と軌道
エネルギー {ε i } を用いると、二次の摂動エネルギーE(2)は次の式で与えられる。

( occ. occ. + vir. 1
221
φi Vpc φ a = ** ∑ +∑ --∑φi Vpc φ a = Eπ(2) + Eσ(2)) i∈π i∈σ , a ε i − ε aa εi − εaocc. vir.

E (2) = ∑∑
i

	
  	
  	
  (1)

ここで occ.と vir.はそれぞれ被占及び仮想軌道を表す。式(1)で定義された E(2)の分解に対応して、

δµは π 及びσ電子の揺らぎの寄与 δµ π , δµσ に分解できる。まず、電子揺らぎの程度を記述するカップリングパラメータλに依存して 0 から E(2)まで変化する溶質­溶媒間相互作用ポテンシャル Eλ(2) を導入する。すると、自由エネルギーδµは Kirkwood の charging 式1δµ = ∫ 0 d λ ∫ d Xρ ( X; λ )d (2) $ &E % X'dλ λ

1= ∫ 0 d λ ∫ d Xρ ( X; λ )1d (2) $ &d (2) $ &Eπ ,λ %X' + ∫ 0 d λ ∫ d Xρ ( X; λ )E %X' = δµ π + δµσ 	
  	
  	
  	
  	
  (2)
dλ
d λ σ ,λで表される。ここで X は全溶媒分子の配置であり、 ρ ( X; λ ) はカップリングパラメータλにおける溶媒の構造 X の出現確率を表す。式(2)を出発点とすれば、通常の ER 法と同様の枠組みで自由エ
ネルギー汎関数を構築することが出来る。従って、 Eπ(2) をエネルギー座標ηにとり、溶液系及び参照系においてηの分布関数 ρπ (η ) , ρπ 0 (η ) をそれぞれ構築することによって、π 電子の揺らぎの寄与
δµ π が計算できる。 δµσ についても全く同様である。

【計算条件】 本研究では、溶質を QM 系で表し、溶媒を MM 系で扱う典型的な QM/MM 計算を用いて、ベンゼン分子の水への溶媒和自由エネルギーに対する π 電子及びσ電子の寄与を計算した。
QM 系の電子状態は実空間グリッドを基底とする Kohn-Sham の DFT により決定し、交換相関エネ
ルギーは BLYP 汎関数を用いて評価した。QM の実空間セルは一辺 9.7Åの立方体とし、各軸方向に 64 個のグリッドを配置した。QM 系の Lennard-Jones 相互作用のパラメータは GROMOS(43A1)
モデルの値を用いた。MM 系として扱う水の溶媒は 496 個の SPC/E モデルの水分子で構成した。
溶媒の配置は分子動力学シミュレーションにより生成し、その時間ステップは 1 fs、溶液の熱力
学条件は 300 K, 1.0g/cm3 に設定した。
【結果と考察】溶液系及び参照系におけるベンゼンに対するエネルギー分布関数 ρπ (η ) , ρσ (η )及び ρπ ,0 (η ) , ρσ ,0 (η ) を図１に示す。ここで、参照系とは先に述べた、溶質の電子密度を固定した中間状態を指す。σ電子と比較して、π 電子は溶媒分子に応答して大きく分極し、安定化に大きな寄与をしていることが分かる。また、溶媒の運動に応答して溶質の電子密度が揺らぐ溶液系において、電子分極が大きいことが分かる。これらの分布関数を用いて求めた自由エネルギー δµ π と δµσを表１に示した。二体的な相互作用の寄与 Δµ は正の値であるが、揺らぎの寄与 δµ による安定化
によって溶媒和自由エネルギーΔµが負の値となることが分かる。また、π 電子は数が少ないにも拘らずσ電子を上回る寄与をしており、ベンゼン分子の水和において、 δµ π が主要な安定化因子となっていることが確かめられた。このように、芳香環を含む系の分子シミュレーショ
ンにおいて π 電子の揺らぎを考慮することは必須である。
現在、ベンゼン環に対する置換基効果を調べるため、いくつかの分子に対して計算を行っている。また、溶液中での瞬間的な電子密度の揺らぎを評価する為に、平均外
場からの外場の変化を摂動とした計算を行っており、これらの結果を発表する予定である。
【参考論文】[1] N. Matsubayasi and M. Nakahara, J. Chem.
Phys. 113, 6070 (2000).
[2] H. Takahashi and N. Matubayasi et al., J.
Chem. Phys. 121, 3989 (2004).
[3] H. Takahashi, et al., J. Chem. Phys. 136,
214503 (2012).

図１ ベンゼンのエネルギー分布関数
表１ ベンゼンの溶媒和自由エネルギー（単位：kcal/mol）ΔµBenzene

1.2

δµ

δµ π

δµσ

-0.9

-0.6

Δµ

Δµ Expt.

-0.3

-0.87

3E08
DFT-MD 自由エネルギー計算によるリチウムイオン電池電解液・添加剤の還元反応解析（物材機構 MANA1, JST さきがけ 2, 京大 ESICB3, 富士フイルム 4）○館山	
  佳尚 1,2,3, 後瀉	  敬介 1,4, 袖山慶太郎 1,3, 奥野幸洋 1,4
DFT-MD free energy calculation study on reduction of electrolytes and additives in Li ion batteries(NIMS WPI MANA, JST PRESTO, Kyoto University ESICB, FUJIFILM Corporation)
Yoshitaka Tateyama, Keisuke Ushirogata, Keitaro Sodeyama, Yukihiro Okuno【序】	  リチウムイオン二次電池では、初回充電時に負極/電解液界面で電解液の還元分解が起こり不動態被膜 (SEI) が形成される。SEI は Li イオン導電性を持つ一方で電子導電性を持たないため電解液の還元分解を抑制し、リチウムイオン二次電池は安定したサイクル特性を示すようになる。エ
チレンカーボネート (EC、図 1) 系の電解液に数 wt%のビニレンカーボネート (VC、図 1) を添加
すると電池特性が向上することが知られており、その役割解明のため EC や VC について分子レベルでの還元分解反応の計算が実施されてきた ([1,2]など)。しかしこれらは、クラスターモデル
による計算で溶媒を顕に取り扱わないものや、溶媒を顕に取り扱ったとしてもエネルギーに基づく議論は行われていないものであった。そこで本研究では、顕に溶媒のダイ
ナミクスを考慮した自由エネルギーに基づく平衡状態および反応自由エネルギー解析を行い、EC や VC の還元後の反応過程を明らかにし、SEI 形成における VC の役割に

図１：VC および EC の分子構造	 ついて知見を得ることを目的とした。
【方法】
	  CPMD コードによる DFT 分子動力学計算を用いた。周期的境界条件を課した一辺 15.24 Å3 の立
方体のスーパーセルを用い、分子数は EC の密度が 1.32 g/cm3 に一致するように調整した。k 点としてはΓ点のみを用い、交換相関汎関数には PBE を用いた。Goedecker タイプ擬ポテンシャルを用い、平面波基底のエネルギーカットオフは 90 Ry に設定した。ダイナミクスについては Nosethermostat を用いて 353 K の NVT アンサンブルを取った。これらの計算条件を用いて、Li+イオン
を含む EC 単独系および VC が添加した EC/VC 系の１電子還元、２電子還元およびラジカル攻撃
反応の平衡構造計算を行った。更に主要な還元分解過程およびラジカル攻撃後の CO2 生成過程については Blue-moon ensemble 法を用いて化学反応自由エネルギー変化の計算も行った。
【結果】
VC および EC の還元分解反応について検討した。EC からは 1 電子還元により開環反応が進行してアニオンラジカル (oE-EC-) を生成し、また 2 電子還元により CO32-、C2H4 や CO が生成した。

一方、１電子還元された VC については、文献[1]で報告されている開環反応よりも、CO 脱離の反応が優位であることが分かった。続いて、還元反応により生成したラジカル種の後続反応としてラジカル種と EC、VC との反応を検討した。oE-EC-は EC とは反応しないが、oE-EC-や VC アニ
オンラジカルは VC と反応し、さらに分解して VC 添加系での特徴である CO2 が脱離することが分かった (図 2)。一方、VC の還元分解により生成したアニオンラジカルは EC や VC との反応が起こらず、反応性に乏しいことが分かった。

図 2.	
 (左)	
 第一原理分子動力学計算とブルームーンアンサンブル法により得られた EC ラジカルアニオンと VC 添加剤の重合過程における化学反応の自由エ
ネルギー曲線と初期、最終および中間状態の構造。この反応が十分起こりやすいことを示している。(右)	
 EC 溶媒中のこの重合反応の最終生成物。これが添加剤導入時の SEI 形成の素材となる。	
 本検討により得られたガス生成物 (CO,






C2H4, CO2) は、EC 系電解液で実験によECり観測されている主なガス生成物 [3]EC-と一致した。また以上の結果から、VCの SEI 形成における本質的な役割は、

VC
SEI

自身が還元分解されることではなく、
図 3.	
 本研究が明らかにした VC	
 添加剤導入による機能の向上した SEI の形成機構。	
 EC ラジカル種を不動態化することで、１電子還元で SEI 素材を形成できる点にあることが分かった[4]（図３）。

[1] Y. Wang, et al. J. Am. Chem. Soc., 124, 4408 (2002). [2] K. Leung, et al. Phys. Chem. Chem. Phys.,12, 6583 (2010). [3] H. Ota, et al. J. Electrochem. Soc. 151, A1659, (2004) [4] K. Ushirogata, K.
Sodeyama, Y. Okuno, Y. Tateyama, J. Am Chem. Soc., in press (2013).

3E09
イオン輸送に伴う液体界面の構造変化の役割(東北大学大学院理学研究科) ○吉川信明，森田 明弘

Structural Fluctuation of Liquid Interface
caused by Ion Transport(Tohoku Univ. Graduate School of Science) ○Nobuaki Kikkawa, Akihiro Morita【序】 親水性の物質が液体界面を通過する際，water finger と呼ばれる水の柱が形成する。
water finger の形成は物質の界面通過速度の決定因子であることが疑われ，MD による water
finger の発見以来[1]，理論的[5, 6]，計算科学的[1-4] 解析が試みられてきた。
我々は近年，イオン輸送を促進する相間移動触媒の MD シミュレーションにより，触媒として働く tetrabutylammonium cation (TBA+) が water finger の切断を促進し，Cl－の輸送自由エ
ネルギーを下げることを見出した[3]。 この結果は water finger の形成切断と輸送自由エネル
ギー曲線との間の関連を示唆する。 そこで我々はさらなる解析のため，water finger の構造を座標化し，その座標を用いた自由エネルギー曲面を用いることを考えた。
【方法】 本解析を行う上で鍵を握るのは water finger の形成切断をどのように射影するかである。 特にこの形成切断は遅いプロセスであることから，拡張アンサンブル法による制御が可能な座標でなければならない。 具体的には，(a) 瞬間の値をトラジェクトリーから厳密に定義できること，(b) 微分（ヤコビアン）の計算ができること，が必要となる。また，(c) 水
和核との区別ができること も重要である。
我々は上記の条件を満たす座標として，water finger
を作る水分子間（または水分子-輸送分子間）の界面
垂直方向の間隔の最大値 zwf を利用する（図 1）
。 zwf

z

max

zwf

は小さいときは water finger が形成している状態，大きいときは切断された状態を射影している。
自由エネルギーの計算法としては，ハミルトニア
ンレプリカ交換 MD 法[7]を用いた。 バイアスポテン
シャルは
2
2
0
�𝑧𝑧wf − 𝑧𝑧wf,𝑖𝑖�
𝑘𝑘B 𝑇𝑇 �𝑧𝑧 − 𝑧𝑧𝑖𝑖0 �
bias
𝑈𝑈𝑖𝑖 (𝑧𝑧, 𝑧𝑧wf ) =
+
�
�
2
2
𝜎𝜎𝑧𝑧,𝑖𝑖
𝜎𝜎𝑧𝑧2wf,𝑖𝑖

を用い，WHAM[8]により界面垂直方向 z と water finger

図 1 water finger を射影した座
標. 赤の分子は水分子，黄色の分
子は輸送分子.

の座標 zwf の 2 次元の自由エネルギー曲面を作成した。
0, 𝜎𝜎𝑧𝑧,𝑖𝑖 , 𝜎𝜎𝑧𝑧wf,𝑖𝑖 はバイアスのパラメータである。 また，𝑧𝑧 ≈なお，i はバイアスの番号，𝑧𝑧𝑖𝑖0 , 𝑧𝑧wf,𝑖𝑖

0 を界面の位置にとる。

【結果と考察】 水気液界面における I－の輸送を解析した例を図 xxx に示す。結果からイオンが水中から気相中へ移動する際，water finger がすぐ切断されるパスと切断されないパスの 2
つが存在することがわかる。両者のパスの自由エネルギーが同程度となる𝑧𝑧 ≈ 20Åでは，両者の間におよそ 4 kacl/mol の活性化障壁が存在する。 この障壁は z のみの自由エネルギー曲

線では観測することができない。
過去に行われた輸送自由エネルギー計算により，water finger が形成切断する位置付近，物
質輸送の向きによるヒステリシスの存在すること[1-3]や摩擦係数が大幅に増大すること[4]等が知られている。 これらの現象は，イオンの移動に 2 種類のパスが存在することや，water fingerの形成切断方向の障壁が z への射影で見えなくなり，代わりに摩擦の項として現れることを考えるとうまく説明できる。
当日は可能ならば液液界面の系について議論する。
【謝辞】 本研究は日本学術振興会特別研究員奨励費の支援を受けて行われた。

[1] I. Benjamin, Science 261 (1993) 1558
[2] K. J. Schweighofer, I. Benjamin, J. Phys. Chem. 99 (1995) 9974
[3] N. Kikkawa, T. Ishiyama, A. Morita, Chem. Phys. Lett. 534 (2012) 19[4] A. Gupta, et al., Phys. Rev. E 78 (2008) 041605[5] R. A. Marcus, J. Chem. Phys. 113 (2000) 22[6] A. A. Kornyshev, M. Urbakh, et al., J. Chem. Phys. 117 (2002) 8[7] H. Fukunishi, O. Watanabe, S. Takada, J. Chem. Phys. 116 (2002) 9058
[8] S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman, J. M. Rosenberg, J. Comput. Chem. 13(1992) 1011
図 2 自由エネルギー曲面.

𝑧𝑧 ≈ 0 が界面であり，z < 0 の領域ではイオンは水相中に存在zwf Åする。実践は山線，点線は谷線を
意味する.

zÅ

3E10
界面における水の核の量子効果（マックスプランク・ポリマー研究所）○永田 勇樹

Nuclear Quantum Effects at the Water-Vapor Interface(Max-Planck Institute for Polymer Research) Yuki Nagata

Using combined theoretical and experimental approaches, we demonstrate that the bondorientation of water at the water-vapor interface depends markedly on the water isotope(H/D) composition. While the interfacial water structures of H2O and D2O are
indistinguishable, the intramolecular symmetry breaking in HDO is directly reflected at thesurface: the OD bonds preferably orient down towards the bulk water, whereas the OH bondtends to orient up into the vapor phase. Path integral molecular dynamics simulations show
good agreement with surface-specific sum-frequency generation (SFG) spectroscopy results,
revealing that the distinct interfacial bond orientations originate from nuclear quantum
effects. The enhanced localization of the heavier D atom leads to stronger hydrogen bonds,
giving rise to OD bonds of HOD molecules preferentially pointing down into the bulk.

Fig. 1.(Left) Pictorial representation of
the angle θ between the OH (OD)bond and the surface normal. (Right)
Axial distribution of average angles<cosθ>. Lines are to guide the eye.
The water molecules illustrate theaverage orientations of HDO (D dark
blue, H light blue) in the different
regions. The positive and negative zaxis points up to the vapor region andto the bulk water region, respectively.

Fig. 2. (a, b) Experimentally measured SFG spectra (red and green) and their fitted curves(blue), (c, d) reconstructed from the fit parameters, and (e, f) simulated. Top panels (a), (c),(e) and down panels (b), (d), (f) display the spectra in the OD and OH stretching regions,
respectively. Filled areas in (c), (d), (e), and (f) represent the free OH or OD bond peak areas.

Reference: Y. Nagata, R. E. Pool, E. H. G. Backus, and M. Bonn, Phys. Rev. Lett., 109,22610 (2012).