4E12

キュミュラント動力学法によるクラスター構造転移の変分的自由エネルギー解析（阪大院基礎工）○重田 育照，桑原 弘幸，乾 智也，

	
  	
 

	
  	
 

馬場 剛史，奥野 克樹，岸 亮平，中野 雅由

Variational free energy analysis on cluster structural transition by cumulant dynamics(Grad. School of Engineering Science, Osaka Univ.) Yasuteru Shigeta, Hiroyuki Kuwabara,
Tomoya Inui, Takeshi Baba, Katsuki Okuno, Ryohei Kishi, Masayoshi Nakano
	 【序】我々はこれまで、量子揺らぎの変数を導入した拡張動力学である量子論的キュ
ミュラント動力学理論を構築し、水素移動反応等の弱い量子効果の解析を行ってきた[1-2]。この方法の本質は、量子論的な期待値であるキュミュラント変数で特徴づけられる分布の動力学を扱う点にある。前回は古典統計力学の問題を取り扱うため、虚時間発展に対する運動方程式の方法、ならびに、実時間運動方程式に温度（運動エネル
ギー）一定の条件から導かれる拘束力である摩擦項を導入する方法で、平衡状態の構
造平均や密度分布を計算する手法を提案した[3]。本研究では、自由エネルギー汎関数の直接的な変分から自由エネルギーを直接決定する手法を提案し、モースポテンシャ
ルで相互作用するクラスターの自由エネルギー面解析に応用する。
【理論】１次元系に対する２次のキュミュラント動力学理論で現れる局所密度は、" $ ( q ! q0 )2 ! 2% ( p ! p0 ) ( q ! q0 ) + # ( p ! p0 )2 %! ( q, p) =exp $ !
'
2h #$ ! % 2
2" #$ ! % 2
#
&
1

(

)

のように ( q0 , p0 , !, ", # ) の組を指定する事で一義的に定義される。プランク定数 h は統
計アンサンブルに関する規格化条件を満たす様に導入されている。ポテンシャル V ( q )に対する自由エネルギーは変数が一つ減り、 ( q0 , p0 , !, " ) の関数として与えられる。

" ( q ! q0 )2 % k T
1
k BT
2
F=
p0 + ! h !
(1 + ln (2!" )) + ( dsV ( s) exp $ ! 2h! ' ! B2 (1 + ln (2!" ))
2m
2
#
&

(

)

これらの変数に関する自由エネルギーは、変分条件 !F !X = 0( X = q0, p0, !, " ) により求められ、例えば調和振動子の問題に対して厳密解を与える。この自由エネルギーは、
運動エネルギー、ポテンシャルエネルギーに関して完全に分離されており、前者の寄与は常に Maxwell-Boltzmann 分布を与えるが、後者の寄与は V ( q ) に依存して変化する。
	  この手法は、３次元 N 粒子系の系にも容易に拡張出来る。今、２体力で相互作用する N 粒子系のハミルトニアンを(

pi2+ ! V ri " rji=1 2mi
i> j
N

H =!

)

とすると、自由エネルギー汎関数は、

p 20i + hTr [! i ]1 N
F=!
+ ! " dsi ds j ! ( si , q 0i , !i ) ! s j , q 0 j , ! j V sij # ! ! i "i + C
2mi
2 i=1
i=1
i> j

(

N

) ( )

で与えられる。ここで、 q i0 , pi0 は座標と運動量の３次元ベクトル、 !i , ! i は座標と運動量の２次のキュミュラントの３! ( si , q 0i , !i ) = ( 2! )

!3/2

!i

!1/2

３行列、C は定数である。また、	
 

T
" 1
%
exp $ ! ( si ! qi0 ) !i!1 ( si ! qi0 )' 	
 # 2&は１粒子密度である。自由エネルギー曲面は、拘束条件付き最小化によって計算される。例えば、距離の拘束条件をつけた自由エネルギーは、

(

)

kij
2
rij ! rij,0pair(i, j ) 2

F' = F + "で与えられる。ここで、 kij はバネ定数であり、十分大きな値の時 rij ! rij,0 となる。
【計算結果・考察】上記の理論を用いて、古典的なモースクラスターの平衡状態での構造解析を行った。２粒子間に働くポテンシャルはモースポテンシャルで記述される

( )

( (

))

((

として、 VM rij = exp !2 rij ! Re ! exp ! rij ! Re)) を用い、 R = 3 とする。このパラメeータでは、７粒子のモースクラスターの場合、２つの安定構造（PBP、COCT）が存在する。これらの安定構造の対称性は PBP の場合は D5h、COCT の場合は C3v である。
初期構造はランダムから出発し、ある程度計算を進めると、安定構造である PBP もしfree energy (/unit)

くは準安定構造 COCT が得られる。
	  図１に、質量 m=1000、温度 kBT=0.01 での PBP 構造（左上 or 右下）と COCT 構造（右上）をつなぐ自由エネルギー面を示した。ここで、反応座標として PBP 初期安定化構造において最近接あるいは最遠方の距離 r13 と r17 を採用した。２つの構造の自図１ 自由エネルギー面

由エネルギー差は 0.276(/unit)であり、PBP から遷移状態の自由エネルギーバリアは約0.41 (/unit)と見積もられた。分子動力学法で自由エネルギー面を求める方法の一つとして熱力学積分法などがあるが、離散化したパラメータに関して平衡状態になるまで計算し、値が収束するまで繰り返さなければならない。本手法では、平衡状態分布の最適化に反応座標１組あたり数十ステップ程度で到達することが可能であり、平衡状態の問題に対して計算コストの面からみて特に優れている手法であるといえる。また、
分子動力学に代表される粒子描像と積分方程式理論に代表される密度描像の利点を兼ね備えた方法であり、多自由度系の新しい概念構築に寄与しうるものと期待される。
References [1] Y. Shigeta, H. Miyachi, K. Hirao, JCP 125, 244102 (2006). [2] Y. Shigeta, T. Matsui, H.
Miyachi, K. Hirao, BCSJ 81, 1230 (2008). [3] 重田他, 第５回分子科学会（札幌）.

4E13
量子古典混合系近似に基づく分子動力学シミュレーションによる溶液中分子内プロトン移動反応の解析：Marcus 理論および実験との対応(名大院工) 小嶋秀和、山田篤志、岡崎 進

A Mixed Quantum-Classical Molecular Dynamics Study of Intramolecular Proton
Transfer Reaction in Solvents(Nagoya University) H. Kojima, A. Yamada and S. Okazaki【緒言】
プロトン移動は凝縮系における基幹的な反応であり、またプロトンの質量が軽いために量子効果が重要である。本研究グループはこれまでに、量子古典混合系近似の枠組みにおいて化学反応の運動方程式の提図 1. MA の分子内プロトン移動反応

案を行い[1, 2]、この手法を H2O 溶媒中のマロンアルデヒド(MA, 図 1)の分子内プロトン移
動反応へと適用してきた[3]。本研究では溶媒種による溶質-溶媒間相互作用の違いによる化学反応機構の相違を検討する事を目的として、極性溶媒である H2O 溶媒中と無極性溶媒である
Ne 溶媒中での両シミュレーションの結果から、反応ダイナミクスに対する量子効果および溶媒効果の検討をするとともに、本計算と同類の反応の実験結果[4]および反応速度理論との比較を行い、その分子論的機構について解析を行った。

【計算方法】
周期境界条件下の一辺 28 Åのセル中に MA 1 個と溶媒分子 721 個を配置した系を H2O
溶媒および Ne 溶媒に対してそれぞれ用意した。MA の力場には八木らの分子軌道計算(MP2/6-31G**)のデータ[5]を用いて IRC に沿った反応経路に零点振動エネルギー補正を加えた有効１次元ポテンシャルのモデルを採用した。分子間相互作用の LJ パラメータには
OPLS-aa を、クーロン点電荷には MP2/6-31G**レベルの ab initio MO 計算による ESP 電
荷の値を用い、カットオフ半径を 13Åとし、12-13Åの区間で switching 関数を用いた。量子古典混合系近似に基づく分子動力学シミュレーション(量子 MD)において量子系の MA の
時間発展は反応物と生成物に分けられた化学反応の描像による時間依存のシュレディンガー方程式を古典系の Ne 溶媒の Newton 方程式と連立して数値的に解く[2]ことで得た。量子系の状態は Surface Hopping 近似で表わし、状態間遷移は fewest switches アルゴリズムで決定した。量子 MD に加えて全古典近似のシミュレーション(古典 MD)を実行し各 5 ns の軌跡
を得た。

【結果】
軌跡データの統計的な解析を行った結果を表 1 に示す。古典 MD では Ne 溶媒中よりも H2O
溶媒中の方で反応速度が大きかったが、逆に量子 MD では Ne 溶媒中の方で反応速度が大きかった。これは相互作用の大きな極性溶媒中よりも相互作用の小さな無極性溶媒中の方が量子効果により反応速度が大きくなることを示している。この傾向は、反応の自由エネルギー差が小さい場合に溶媒の比誘電率が大きくなるに従い反応速度が小さくなるという溶液内プ
ロトン移動反応の実験結果の傾向[4] と一致している。また古典 MD の結果は遷移状態理論
による反応速度の予測と、そして量子 MD の結果は Marcus 理論による予測と定性的に一致した。
これらの傾向の分子論的機構を得るために解析を行った。量子 MD の反応機構でトンネル
移動反応は両溶媒において大きな割合を占めており、この機構の寄与により古典 MD と量子
MD とで二つの溶媒中での反応速度の逆転が起こっている。これは溶媒中の反応ポテンシャ
ル曲線の形および変化から説明できる。H2O 溶媒中では極性分子の溶媒和により反応物側のポテンシャル井戸が大きく安定化されて非対称型のポテンシャル曲線になり、反応の機会は
大きな熱揺らぎにより反応物側が不安定化した時である。これに対し、Ne 溶媒中では溶媒和による反応物側の安定化は殆どなく、反応ポテンシャル曲線は反応物側と生成物側の安定性が同程度で対称的であり、状態間の振動エネルギー準位が近く常にトンネル移動反応が起こり易い状態にある事と、熱揺らぎにより低いバリアが出現する事で波動関数の重なりが大きくなる事が、トンネル移動反応を促進させる。古典力学的にも起こりう
る熱活性反応とバリア消失時の反応をみると、これらの和は、量子 MD

表 1. Ne および H2O 溶媒中の 5 ns 分の反応とその反応機構では両溶媒で零点エネルギーによる反応の促進が見られたが、両 MD に
おいて H2O 溶媒中の方で反応回数が大きかった。つまりはこれらの反応機構では溶質-溶媒間相互作用による反応促進が支配的であった。

【参考文献】[1] A. Yamada and S. Okazaki, J. Chem. Phys. 124, 094110 (2006). [2] A. Yamada and S.
Okazaki, J. Chem. Phys. 128, 044507 (2008). [3] A. Yamada and S. Okazaki,in preparation. [4] K. S. Peters
and G. Kim, J. Phys. Chem. A 108, 2598 (2004). [5] K. Yagi et al., J. Chem. Phys. 115, 10647 (2001).

4E14

水の分子間振動による電子の挙動と赤外強度(静岡大教育)鳥居肇

Behavior of Electrons Induced by IntermolecularVibrations in Water and Infrared Intensities(Shizuoka University) Hajime Torii
序
液体の水の赤外・遠赤外（テラヘルツ）スペクトルには，幾つかの特徴が知られており，そのうちの１つが，200 cm-1 付近のバンドの存在である。これは，水素結合 (O..H) 伸縮（分子の自由度としては並進）に由来するが，電気的に中性な分子の並進は分子間相互作用無しには赤外強度を持ち得ないので，このバンドの赤外強度は専ら分子間相互作用に由来することになる。最近，このバンドの赤外強度が，分子並進に伴う分子間の電子移動（電荷フラックス）に由来することを，水 28 乃至 30 分子から成るクラスターを対象とした電子密度微分の解析と，その結果を用いた液体水のスペクトルのシミュレーションから明らかにした [1]。
このスペクトルシミュレーションでは，分子並進 (Xζ) および回転 (Ωζ) に関わる双極子微分として，クラスター内で四面体的な水素結合を形成している水分子（計 16 分子）についての計算値の平均［∂μx/∂Xx = 1.91 D Å-1, ∂μy/∂Xy = –1.76 D Å-1, ∂μy/∂Ωx = –∂μx/∂Ωy = 1.57 D (rad-1)］
を用いており，200 cm-1 付近およびそこから~1000 cm-1 までの領域のスペクトルの計算結果は，
実測スペクトル [2] の特徴を良く再現している。しかし，それより低振動数の領域のスペク
トル強度に対する理解を得るためには，液体中に存在する水素結合の欠陥（四面体構造からのずれ）など構造の揺らぎの影響を考慮する必要があると考えられる。また，上述した双極子微分のうち∂μ/∂Ωx および∂μ/∂Ωy は，dynamical dipole moment とも呼ばれるものであり [3]，孤立分子では永久双極子モーメントと一致するが，水素結合に与る水分子については，上述
したように（双極子モーメントが増大するにも関わらず）むしろ孤立分子より小さい値が得られている。1000 cm-1 付近までの領域におけるバンドの相対強度を考える上では，この∂μ/∂Ωxおよび∂μ/∂Ωy の減少のメカニズムも明らかにする必要がある。本研究では，その２点について，理論的検討を行った。
分子回転に関わる双極子微分に対する水素結合形成の影響
水 30 分子から成るクラスターの中心付近に位置する分子について，y 軸まわりの分子回転に関わる電子密度微分の水素結合形成による変化δ(∂ρ(el)(r)/∂Ωy) [≡ (∂ρ(el)(r)/∂Ωy)cluster – (∂ρ(el)(r)/∂Ωy)isolated] を B3LYP/6-31+G(2df,p)レベルで計算した結果を図１に示す。１次元の累積積分プ
ロットから，電子密度が図内左寄りに移動していることがわかる。この電子密度の移動に起因する双極子微分は∂μx/∂Ωy = 0.181 D rad-1 であり，孤立分子の∂μx/∂Ωy = –1.890 D rad-1 を部分的に打ち消して大きさの減少に寄与している。ここで，中央分子の重心を中心とした半径 1.4
Å の球面で空間領域を二分すると，双極子微分に対する球内の電子密度の寄与は∂μx/∂Ωy =–0.042 D rad-1 でむしろ負符号であり，したがって打ち消しの効果は球外の電子密度に由来すると結論づけられる。２次元

-2

-2

6

× 10 a0 rad

× 10 a0 rad

プロットから，中央分子の水素原子の運2動（反時計回り）によって，水素結合受1

容体分子の酸素原子の電子が引きずられ

-1

0
-2

て移動していることがわかる。同様のことは，x 軸まわりの分子回転に関わる電子3

密度微分 δ(∂ρ(el)(r)/∂Ωx)でも見られ，これ-3素結合構造揺らぎの影響

-2
-2

-1

× 10 rad

-1

-1

0
-6

が打ち消し効果の主要因と考えられる。
分子並進に関わる双極子微分に対する水-1

図１：水 30 量体クラスターの中心付近に位置する分子の電子密度微分 δ(∂ρ(el)(r)/∂Ωy) の２次元マップ（等高線図）および１次元プロット（黒実線）。

分子並進による分子間の電子移動は，

青点線は１次元プロットの累積積分。図の横軸が x水素結合を介して引き起こされるため，であり，中央の分子が図平面内で反時計回りに回転

水素結合長に依存すると考えられる。そしている。

こで，水素結合方向の分子並進に関わる

長 1.8 Å 程度で ≅ 0.1 e Å-1 であり，水素結合距離が長くなると減少することがわかる。回帰直線から求められる水素結合長の閾値は 2.37 Å であり，液体中に存在す
る水素結合の欠陥との間で，パラメータ12-2を図２に示す。|∂qJ/∂X| の値は，水素結合|∂qJ /∂X| / 10 e Å

スター内の分子を対象に計算した。結果-1

分子部分電荷の微分 (∂qJ/∂X) を，水クラ

108ー値のスムーズな接続が可能となる。

1.7

この結果を用いて 液体水のスペクトル

1.8

1.9

Hydrogen bond length / Åのシミュレーションを行ったところ，前

図２：水 28 量体および 30 量体クラスター内の，四報の場合 [1] と同様，~200 cm-1 のバンド

面体的な水素結合を形成している水分子を対象にの存在を再現することができた。また，計算した，水素結合方向の分子並進に関わる分子部

それより低振動数の領域のスペクトル強度に対する分子並進と回転の寄与の程度
は，前報の場合と異なっており，その振分電荷の微分。赤: 水 28 量体，青: 水 30 量体，○：
水素結合供与分子の並進，●：水素結合受容分子の並進。

動数領域においては水素結合欠陥の影響が大きく見られることが示唆された。
[1] H. Torii, J. Phys. Chem. B 115, 6636 (2011). [2] J. E. Bertie and Z. Lan, Appl. Spectrosc. 50, 1047
(1996). [3] A. Pasquarello and R. Resta, Phys. Rev. B 68, 174302 (2003).

4E15

トランスアゾベンゼン n*/π* 励起による光異性化反応機構に関する理論的研究(1 北大院総合化学、2 北大院理) 原渕祐 1・関奈々美 1・石井萌 2・
中山哲 2・野呂武司 2・武次徹也 2
Theoretical study on photoisomerization mechanism of trans-azobenzene(Hokkaido Univ.) Yu Harabuchi, Nanami Seki, Moe Ishii, Akira Nakayama, Takeshi Noro,
Tetsuya Taketsugu【緒言】
アゾベンゼンはシス-トランス光異性化を示す典型的な分子であり、その光化学特性がナノスケール光デバイスとして利用可能であることから、基礎的観点からだけではなく工業的応用においても大きな可能性を有し、理論・実験双方面
より多くの研究が行われてい図１：アゾベンゼンのトランス－シス光異性化経路

る。初期の実験研究からは、
S0→S1(n*)励起では NNC 結合角が反転することによる inversion、S0→S2(*)励起では NN 結
合周りのフェニル基の回転による rotation と異なる機構で反応が進行すると考えられてきた(図１)[1]。Fujino らは、時間分解ラマン分光実験[2]と過渡吸収スペクトル[3]に基づき、トラ
ンスアゾベンゼンは S2(*)励起後~110 fs という超高速で平面構造の S1(n*)に失活し、S2(*)励起と S1(n*)励起ではどちらも同様の経路を辿ることを示唆した。また、基底状態と S1(n*)状態において NN 伸縮の振動数がほとんど変化しないことから、トランスアゾベンゼンの光異性化反応では、NN の二重結合性が保持され平面構造が保たれると報告している[2]。一方理論研究では、反応に関わる励起状態のポテンシャル曲面の解析[4]においても ab initio 分子動力学計算[5]においても、トランスアゾベンゼンの光異性化反応では n*励起*励起ともに平面構造が崩れて反応が進行する rotation の経路が有利であると報告されている。また、*励起後の光異性化過程では、2 電子励起状態2*2 が重要な役割を担うと報告されている[6]。
本研究では、トランスアゾベンゼンに対し、電子基底状態および S1(n*)状態における構造最適化と振動解析を行い、実験により測定された NN 伸縮の振動数と光異性化反応経路を理論計算の立場から吟味することを目的とした。また、S2(*)励起後の反応機構については、CASPT2 レベルの電子状態計算により最急降下経路を計算することで、これを議論した。

【計算方法】
構造最適化と調和振動解析の ab initio 計算には MP2, QCISD(T), CASSCF, CASPT2 法を適用し、DFT/TDDFT 計算では、B3LYP, CAMB3LYP 法を用いた。基底関数には Sapporo-DZP、Sapporo-DZP+diffuse を用いた。また、ポテンシャル曲面の非調和性を考慮した振動数を見積

もるため、B3LYP(CAMB3LYP)/6-31G*計算に基づき４次力場(QFF)を生成し、cc-VSCF 計算を行った。S2(*)励起後の最急降下経路計算は、CASSCF(6,4)/DZP を参照関数とする CASPT2法を用い、
*状態から n*状態への超高速失活過程について議論した。DFT 計算と CASSCF

計算には GAMESS を用い、CASPT2 と QCISD(T)計算には MOLPRO2010 を用いた。
【計算結果】
全ての計算手法で、n*状態の平面極小構造における振動数は、NN 伸縮の振動数が基底状態平衡構造よりも高くなるという結果が得られた。n*励起に伴い N 原子の孤立電子対から電子が励起するため、CNN 結合角が広がり、更に CN 結合が短くなる一方、NN 結合長は
ほとんど変化しない。図２に CASPT2 計算の結果を示すが、この構造変化に伴い、CN 対称
伸縮振動の振動数がいったん増加し、同じ対称性の NN 伸縮振動モードと混合して NN 対称
伸縮振動の振動数を増加させることが明らかとなった。
当日は、振動解析の結果と*励起状態からの最急降下経路計算の結果の詳細を報告し、それらの結果を踏まえたうえでトランスアゾベンゼンの光異性化反応機構について議論する。

図２：基底状態・励起状態における NN 伸縮・CN 伸縮の振動数変化【参考文献】[1] H. Rau and E. Luddeke, J. Am. Chem. Soc., 104, 1616 (1982).
[2] T. Fujino and T. Tahara, J. Phys. Chem. A, 104, 4203 (2000).
[3] T. Fujino, S. Y. Arzhantsev, and T. Tahara, J. Phys. Chem. A, 105, 8123 (2001).
[4] T. Ishikawa, T. Noro, and T. Shoda, J. Chem. Phys., 115, 7503 (2001).
[5] Y. Ootani, K. Satoh, A. Nakayama, T. Noro, and T. Taketsugu, J. Chem. Phys., 131, 194306 (2009).
[6] I. Conti, M. Garavelli and G. Orlandi, J. Am. Chem. Soc., 130, 5216 (2008).

4E16
蛍光性分子の光化学動力学（1 上智大院理工, 2 レブデフ物理研 3 京大福研セ）村上龍大 1，Kondorskiy Alexey2，石田俊正 3，南部伸孝 1

Photochemical dynamics of fluorescent molecule（1Sophia Univ., 2Lebdev Phys. Inst., 3Kyoto Univ.）
Tatsuhiro Murakami1， Alexey Kondorskiy 2，Toshimasa Ishida3，Shinkoh Nanbu1【序】高選択的な蛍光を示す化合物は、化学、環境学、生物科学の分野で様々なアプリケーシ
ョンとして興味が持たれている。高感度で選択的な蛍光プローブの理論的設計を行うためには光吸収、発光、波長、ストークスシフト、量子収率、エネルギー緩和、無輻射過程などの特徴的な情報を持つ光物理化学過程の解明が必要となる。しかしながら、いまのところ蛍光の量子収率をコントロールする光物理化学過程は十分に理解されていない。そこで本研究として、蛍光性分子について量子効果を取り入れた励起状態 ab initio MD 計算を行い、光物理化学課程を探索した。
【手法】
朱-中村公式に基づく trajectory surface
hopping 法(ZN-TSH 法)により励起状態 ab initio

ab initio MD の方法

MD 計算を行った。励起状態 ab initio MD 法は、初期構造（基底状態からの垂直励起）

電子状態を CASPT2/CASCSF 法により量子力学に基づき計算、原子核に働く力を求め、速度ベルレ法を用い、古典力学に基づき原子核を運動させていく“on-the-fly”の方法である。初期座標と初期運動量は、ゼロ点振動エネルギーの振動

原子核に働く力の計算(ab initio 計算による電子状態計算）

波動関数の形であるウィグナー分布を基に決定繰り返しした。ゼロ点振動波動関数は基底状態において
CASSCF 法によって得られた解析的ヘッセ行列

原子核の運動の計算を用いて計算した。励起状態から基底状態への(古典力学）

遷移が予想されるため、ZN-TSH 法を用い毎時
間ステップに非断熱遷移が起こるかどうかをチ
ェックしながら、遷移が起こる際には、ZN の非断熱遷移のチェック
図 1. ab initio MD 法の概略模式図

公式を用いて非断熱遷移を正確に取り扱う[図 1]。ZN-TSH 法は、非断熱遷移を起こす可能性のある領域に trajectory が来たときに、散乱理論の式を使用し遷移確率を求めるという手法であり、TSH
法の散乱問題において、半古典 ZN 理論と組み合わせることにより、古典的禁制遷移も可能にすることができる。また、ZN-TSH 法では、有名な fewest switches TSH 法で発生する overcoherence
問題も生じない。

さらに、ZN-TSH 法では欠落している位相情報を Herman と Kluk により提唱された凍結
ガウス関数の(以下 HK-FG)時間発展法を応用し、量子効果を取り入れた。Herman と Kluk の方法
は、核の波動関数をガウス関数で展開した波束を古典力学に従って伝播させるという方法であり、個々のガウス関数は時間発展の最中で変化しないので凍結ガウス関数と呼ばれる。核の波動関数は次式で表される。

Φ m (t )i

dq 0 dp 0

exp  Sg (q , p ) C∫ ( 2p  ) = g (q , p ) Φ (t 0)


N

t

t

q0p0 ,tq0p0 ,t

0

0

(1)

m

時刻 t における全波動関数は、個々のガウス関数の時間発展により得られた波動関数 Φ m (t ) より
以下の通り求められる。

Ψ (q, r; t ) =∑ Φ n (q;t )φn (q, r)n

そして、Tatchen らが最近提案した“poor person’s” HK-FG(PP-HK-FG)時間発展法を利用することで、式(1)の係数 C の部分で表れるヘッセ行列を求める必要がなくなる。ZN-TSH 法では、式(1)で表れる初期座標・初期運動量、毎時間ステップの座標・運動量・運動エネルギー・ポテンシャルエネ
ルギーを求めているので、この手法と組み合わせることにより、自己相関関数 A(t ) = Ψ (t =0) Ψ (t )を求め、この自己相関関数をエネルギーと時間の間で Fourier 変換することにより理論吸収スペク
トルを設計した。
【結果と考察】
蛍光性分子として、発光分析試薬や赤色 LED に興味が持たれているインドリルマレイミ
ド(IM)[図 2]を主に取り上げた。IM は極性非プロトン性溶媒のアセトニトリル中で著しいストー
クスシフトを伴う蛍光を示す分子である。我々は最近、ZN-TSH 法により、IM 分子のインドール
基とマレイミド基から脱プロトン化した 2 価のアニオン種(IM2-)が発光化学種であると帰属した 1。
IM2-の S0-S1 垂直励起エネルギーは 37840 cm-1(CASSCF/cc-pVDZ)であり、ZN-TSH 計算+PP-HK-FG
時間発展計算によって得られた理論吸収スペクトルとのピーク位置がほぼ一致した[図 3]。詳細は公演にて発表する。

0.05

O

N
H
図 2. IM の分子構造

IntensityO

H
N

37919 cm

0.04

0.03

0.02

0.01

0.00

30000

40000

50000

Wavenumber/cm【参考文献】

-1

60000

-1

図 3. IM2-の理論吸収スペクトル[1]. T. Murakami et al., Phys. Chem. Chem. Phys., 14, 11546 (2012).

4E17
完全固体・液体 v3(法政大生命科学) 片岡 洋右，山田祐理

Thermodynamics of 3-Phase Equilibrium in Argon
Based on a Perfect Solid and a Perfect Liquid (v3)
(Hosei University) Yosuke KATAOKA and Yuri YAMADA
1 初めにアルゴンにおける三相の間の相転移を議論できる固体と液体を単純化したモデルを与える。
これを完全固体・液体と呼ぶ[1.2]。分子間にレナード・ジョーンズ（LJ）関数を仮定する。

  12   6 u( r )  4      
 r   r ポテンシャルの深さである をエネルギーの単位とし, 直径に相当する を長さの単位とする。
2 基本的考え方
低温で LJ 固体は FCC 格子作る。一定密度で温度を上げ、ポテンシャルエネルギーの平均値 Ue は調和振動子を仮定すると温度に比例する温度効果項を持つ。圧力 p も調和振動子を仮定すると同様に低温での値に、温度効果項を加えて有限の温度における値を表すことができる。
3 状態方程式の構成
LJ 系について分子動力学シミュレーション(MD)を行い、体積 V ごとに、ポテンシャルエネルギ
ーの平均値 Ue と圧力 p の温度変化をプロットすると、図 1,2 のようにほぼ直線が得られる。
-230(N/V)=0.973
Lines: Fitting
Dots: MD

-3(N/V)=0.97
Lines: Fitting
Dots: MD

-3

25
20
Liquidp/(/3)(Ue/N)/

-4
-5
-6

10
5

-7

0

Solid
-8

Liquid
15

0

1

2

3
T/(/k)

4

-5

5

Solid
0

1

図1

MD の Ue と p を T の 1 次式で表した。

UekT U e (0K)
 a(U e )

N

N

2

3

4

p
kT p(0K) a( p)

3
 /

 /3

この係数と定数項を図 3-6 に示した。定数項は数密度の２項の多項式で表した。
6
U e,s (v,0 K )1   12
 1
 6 1  4  12  1   2 (Solid)N 128  v
 5 v

U e, f (v,0 K )
N

5

図2

T/(/k)  18
3 
 1.0692052  1.5 6  8.2  (Liquid)
v
v 


12

2.5
-3

as=1.35*(N/V)/

solida(p)*(V/N)/ 3

a(Ue)

2

Solid

1.5

a(p)*(V/N)/3=1+6(N/V)/-3

10

1
Liquid86liquid4a(p)*(V/N)/3

2=1+4.9(N/V)/-3

0.5
0

-3

af=(N/V)/

00

0.2 0.4 0.6 0.8 1(N/V)/-3
Ep(0K) 16:59:46-2

1.2 1.4 1.6

0.2 0.4 0.6 0.8 1(N/V)/-3

図3

1.2 1.4 1.6

図42012/06/19

15

250
Uef(0K)Uef function of v
Ues(0K)Ues function of v
Ues by sum

ps(0K)
pf(0K)-dUes(0K)/dV-dUef(0K)/dV

200p(0K)/(/3)

10(Ue(0K)/N)/ 

05
Liquid0-5150Liquid

100
50
0

Solid

Solid

-50

-10
0

0.2 0.4 0.6 0.8 1(N/V)/-3

1.2 1.4 1.60

図5

0.2 0.4 0.6 0.8 1(N/V)/-3

1.2 1.4 1.6

図6

4 完全固体・液体の状態方程式

3
NkT  U e (V ,0 K )  g (V ) NkT
2
NkTd dg (V ) p(V , T ) U e (V ,0 K )  f (V ) NkT   NkT ln  kT VdV dV U (V , T ) 

熱力学的状態方程式を満たすためには、圧力の式の最後の補正項を加えておく必要がある。
5 アルゴンの 3 相平衡
ほかの熱力学関数は標準的手順で定義できる。与えられた温度に対しギブズエネルギーの相平衡の条件式を数値的に解けば、相平衡を求めることができる。図 7,8 に相境界の圧力―温度曲線と温度―数密度曲線を示す。今回得られた結果は巨視的実験値や MC,MD の結果に良く対応する。
105
10

5

Solid

4

1000
LiquidT/(/k)

100

p/atm

Lines: EOS v3
Dots: MC & MD

4

10
1

3
2
Gas

0.1
Gas

0.01

1

Lines: EOS v3
Dots: exp, MC & MD

0

0.001
0

100

200

300
T/K

400

500

Liquid

600

図７

Solid

0

0.2

0.4

0.6 0.8(N/V)/-3

1

1.2

文献 [1] Y. Kataoka and Y. Yamada, J. Comput. Chem. Jpn., 10, 98-104 (2011).
[2] Y. Kataoka and Y. Yamada, J. Comput. Chem. Jpn.. in press.

1.4

図8

4E18

Multi-product 展開を用いた Reptation Monte Carlo 法の開発
（横市大院・生命ナノ）北 幸海，立川 仁典

Development of Reptation Monte Carlo method with higher orderpropagator based on Multi-product expansion(Yokohama-city Univeristy) Yukiumi KITA, Masanori TACHIKAWA	 【 緒 言 】 	 
	  量子モンテカルロ（QMC）法[1]は、原子や分子そして固体を問わず適用可能な高精度第一原理法として注目されている。特に、時間依存 Schrödinger 方程式の虚時間発展を利用した拡散モンテカルロ（DMC）法により、従来の量子化学的手法と比較して遥かに高い精度で変分エネルギーの解析が可能である[2]。しかしながら、(i) 定量的に信頼できる結果を得るためには虚時間ステップに関する外挿が必要不可欠である、(ii) 配置数変動による誤差（population
control bias）が不可避である、といった問題も知られており、これらの問題は大規模系や遷移金属等、電子数の多い系では特に深刻である。さらに、(iii) mixed estimator という特殊な期待値の算定方法を用いるため、系のハミルトニアンと非可換な物理量の算定が原理的に難しく、
これは QMC 法の適用範囲を著しく狭める要因の１つとなっている。
	  Reptation Monte Carlo (RMC) 法[3]は、DMC 法と同様、虚時間発展を利用した QMC 法の一つであるが、原理上、配置数変動による誤差を伴わず、ハミルトニアンと非可換な物理量を
厳密な表式で算定可能な QMC 法である。また近年 Chin によって開発された Multi-product（MP）展開法[4]は、指数型演算子に対する高次展開法の１つであり、二次のプロパゲータのみを組み合わせるだけで、偶数次の高次プロパゲータを実装可能な手法である。したがって、
MP 展開法による高次プロパゲータを実装した RMC 法を開発する事で、(i)虚時間ステップの外挿が不要であり、(ii)配置数変動による誤差を原理上伴わず、かつ(iii)ハミルトニアンと非可換な物理量の厳密算定が可能な、新しい量子モンテカルロ法の開発が可能となる。
	  そこで本研究では、Multi-product 展開による高次の虚時間プロパゲータを用いた Reptation
Monte Carlo 法の定式化、およびプログラムへの実装を行い、幾つかの小分子の振動状態を例にその有効性を検証した。
【 方 法 】 	 	
  時間依存 Schrödinger 方程式の虚時間発展を利用した QMC 法では、任意の試行波動関数! T に虚時間プロパゲータを作用させる事で系の正確な固有状態 ! 0 を抽出する：!
$ Ĥ

! 0 = lim e 2 ! T .
! "#

(1)

ここで Ĥ は系のハミルトニアン、! は虚時間を表す。Reptation Monte Carlo 法は、以下で定義される仮想的な分配関数 Z0 を導入し、この分配関数をメトロポリス法によってサンプリング
する事で、系の正確な波動関数とエネルギー固有値を求める方法である：

Z 0 ! " 0 " 0 = lim " T e %! Ĥ " T ,
! #$

(2)

ここで ! 0 は系の正確な固有関数を表す。RMC 法では、経路積分形式で虚時間プロパゲー
タを短時間分解した後、短時間近似によって得られる近似的な分配関数を用いる。例えば、二次の短時間近似を用いた分配関数は以下のように表される："!
N #1# {" ( x j )+" ( x j+1 )}% N #1(Z 0 = , ! , ! T (x0 ) '$ g ( xi , xi+1; "! )* ! T (x N ) + $ e 2
dx0 !dx N .
& i=0
)
j=0

(3)

	
 
ここで N は虚時間プロパゲータの分割数、 !! (= ! / N ) は虚時間ステップ、 ! ! " #1
T Ĥ " T は局
2
所エネルギー、g は ! T に従う分布を生成する短時間プロパゲータである。
	  一方、近年 Chin によって開発された Multi-product 展開法は、以下のように二次のプロパゲ
ータ G2 を組み合わせる事で、偶数次の高次プロパゲータ G2n を表す方法である：" !! %
G2n (!! ) = ( ciG $+ " (!! 2n+1 ),'# ki &
i=1
n

ki
2

ki2
.
ci = # 2
2
j=1("i) ki ! k jn(4)

Zero-point energy [cm-1]

ここで{ki}(i=1~n)は同一の値を含まない限り任意の値である。本研究では MP 展開法を用いて、2~24 次の高次プロパゲータの実装を行った。
	 【 結 果 と 考 察 】 	
 	
  開発した手法の有効性を検証するた
2205
6o=1 x 10-6 [K-1]
め、H2 および H2O 分子の振動状態の解
6o=1 x 10-5 [K-1]2200
析を行った。RMC 計算における配置空6o=3 x 10-5 [K-1]6o=6 x 10-5 [K-1]間の定義には、分子の平衡構造で定義
2195
6o=1 x 10-4 [K-1]
された基準振動座標を用い、ポテンシ
Numerial exact
ャ ル エ ネ ル ギ ー 曲 面 (PES) は 、
2190
Born-Oppenheimer 近似の下で第一原理
2185
計算から得られた grid 点を多次元スプ
ライン補完することで生成した（計算2180
レベルは H2 では Full-CI/aug-cc-pVTZ、
2175
H2O では CCSD(T)/ aug-cc-pVTZ）。
2468
10 12 14 16 18 20 22 24	  Figure 1 に開発した手法による H2 分
Multi-product order
子の ZPE の虚時間ステップ依存性を示 Fig. 1. The multi-product order dependence on the zero-point energy(ZPE) of H2 molecules obtained with reptation Monte Carlo
す。Figure 1 より、最もプロパゲータの calculations. The error bar in each symbol means the standard
精度が低い NMP = 2 の場合には、非常に deviation of ZPE. The numerical exact value obtained with
Störmer-Levy method (black solid line) is also shown.
小さな時間ステップ( !! = 1"10 #6 K #1 )を用いなければ、数値的厳密解と同等の結果は得られないことがわかる。一方、MP 展開による高次プロパゲータを用いることで ZPE は系統的に改善し、NMP ≧22 では使用した全ての時
間ステップで厳密解と同等の結果が得られることがわかった。発表当日は、H2O 分子の振動状態解析の結果についても報告を行う予定である。
【参 考 文 献 】[1] B. L. Hammond, W. A. Lester Jr. and P. J. Reynolds, “Monte Carlo Methods in Ab Initio Quantum
Chemistry” (World Scientiﬁc, 1994). [2] Y. Kita, et al., J. Chem. Phys. 131, 134310 (2009); Y. Kita etal. J. Chem. Phys. in press (2011). [3] S. Baroni and S. Moroni, Phys. Rev. Lett., 82, 4745 (1999). [4]S. A. Chin, e-print arXiv:0809.0914v2 (2009); R. E. Zillich, J. M. Mayrhofer, and S. A. Chin, J. Chem.
Phys., 132, 044103 (2010).

4E19

大規模分子の超並列計算のための MPI/Open-MP ハイブリッド
並列 RI-MP2 アルゴリズム（理研 AICS）河東田 道夫，中嶋	  隆人

MPI/Open-MP hybrid parallel RI-MP2 algorithm for
massively parallel computations of large molecules(RIKEN AICS) Michio Katouda, Takahito Nakajima【序】近年のスーパーコンピュータや PC クラスタでは、マルチコア CPU を複数搭載したノードを高速なインターコネクトで接続したマルチコア超並列クラスタシステム構成が一般的となっている。現
在では京コンピュータのような数万ノード・数十万 CPU コアで構成された 10 ペタフロップス級の演算性能を持つスーパーコンピュータも登場している。これらのシステムを活用し、ナノ分子の機能のデザインや生体系の現象解明のために巨大分子の量子化学計算を行うためには、マルチコア超並列ク
ラスタシステムの性能を十分に引き出せるアルゴリズムとプログラムの開発が重要となっている。特にナノマテリアルや生体分子といった巨大分子内では、ファンデルワールス力などの非結合相互作用が重要な寄与をするので、これらの相互作用を高精度に記述することが可能な Møller–Plesset 2 次摂動(MP2)法に基づいた計算が必要となる。本研究では、巨大分子の RI-MP2 計算を京コンピュータな
どのマルチコア超並列クラスタシステム上で高速・高並列に行うことを目的として、RI-MP2 法の
MPI/Open-MP ハイブリッド並列アルゴリズムの開発を行った。
【方法】発表者らは以前の研究において、大規模分子の MP2 計算を PC クラスタ上で行うことを目的
として、RI-MP2 法の並列計算アルゴリズムおよびプログラムの開発を行った[1,2]。しかしながら、以前に開発した RI-MP2 並列計算アルゴリズムでは、MPI あるいは TCP/IP ソケットを利用したプロ
セス並列化のみを行ったため、マルチコア超並列クラスタシステムのノード内共有メモリ・通信イン
ターフェースを効率よく使って高並列性能を出すことが困難であった。また、プロセス並列のタスク
分割を占有軌道を用いて行なっていたため、並列度が占有軌道数で制限され、プロセス数の増加に対して並列性能が頭打ちになることも問題となっていた。そこで、本研究ではノード間の並列化を MPI、ノード内の並列化を Open-MP で行う MPI/Open-MP ハイブリッド並列化を行い，仮想軌道を並列分割軸に使うことで、これらの問題を解決することが可能な並列アルゴリズムを開発した。
図 1 に本研究で開発した RI-MP2 計算の並列アルゴリズムの概要を示す。この図に示す 3 つのステップの処理について、最外のループを MPI でプロセス並列を行なっている。また、ステップ 2、3 の処
理では、以前に開発したアルゴリズムでは占有軌道のループを MPI 並列化を行なっていたが、本アルゴリズムでは仮想軌道のループが最外になるようにループの順を入れ替えて MPI 並列化を行った。一般に仮想軌道の数は占有軌道の数より多くなるので、以前のアルゴリズムと比べて並列プロセス数を増やした場合でもより高い並列性能を達成することが可能である。MPI 並列化された最外ループの内側で行われる処理については Open-MP でスレッド並列化を行っている。RI-MP2 計算では計算コスト
的にボトルネックとなる 3 中心積分の積分変換や 4 中心積分生成計算は行列-行列積計算の形で処理を行うことが可能である。そこで、本アルゴリズムを実装する際には、Intel Math Kernel Library などの最適化 BLAS ライブラリに実装されているスレッド並列版の DGEMM を用いることで高いスレッド並列性能および演算性能を容易に達成することが可能である。
[Step 1]Loop µ (MPI parallelization)
Evaluation of three-center AO ERIs (µν|l)(OpenMP parallelization)
1/3 and 2/3 index transformation of (µν|l)by threaded DGEMM in BLAS (OpenMP parallelization)(ia l ) = !C! a !Cµi (µ! l )!

[Step 3]Loop batch of virtual orbital A (MPI parallelization)
Read B ia (a ! A)n

Loop bProc
Send B ia to Myrank- bProc nodenReceiveµ

Loop a ! A (MPI parallelization)
Loop b ! Myrank + b (MPI parallelization)

End Loop µ

Proc[Step 2]
Evaluation of two-center AO ERIs (l|m)(OpenMP parallelization)
Cholesky decomposition and inversion of (l|m)
by threaded LAPACK (OpenMP parallelization)
Loop virtual orbital batches A (MPI parallelization)
3/3 index transformation of (ia|l)by threaded DGEMM in BLAS (OpenMP parallelization)

Bnia = " (ia l ) (l n )!1/2

Bnjb from Myrank+ bProc node

Evaluation of four-center MO ERIs by threaded
DGEMM in BLAS (OpenMP parallelization)(ia | jb) = ∑ Bnia BnjbnEvaluation of MP2 correlation energy E(2)(OpenMP parallelization)
End Loop b
End Loop a
End Loop bProc
End Loop A(a # A)l

End Loop A

図 1. MPI/Open-MP2 ハイブリッド並列 RI-MP2 アルゴリズム【結果】今回開発したアルゴリズムを分子科学計算プログラム NTChem に実装し、理研 RICC 超並列
クラスタシステムで並列性能のテストを行った。図 2(a)にタキソールの RI-MP2/cc-pVDZ 計算の CPU
コア数に対する高速化率を示す。MPI/Open-MP ハイブリッド並列化によりフラット MPI 並列化の場
合より高速化率が大きく向上した。また、MPI 並列化の分割軸を仮想軌道にしたことにより、分割軸を占有軌道とした時よりも高速化率が大きく向上した。図 2(b)はバッキーキャッチャーC60@C60H28 のRI-MP2/cc-pVTZ 計算の高速化率を示す。512 ノード 4096CPU コアを用いた場合でも 2621 倍の高速化率を達成した。また、512 ノード 4096CPU コアを用いた場合、C60@C60H28 の計算は 15 分で終了した。本研究で開発したアルゴリズムおよびプログラムを用い、京コンピュータなどの超並列計算機上で数千 CPU コアを用いて計算を行うことにより、300 原子程度のサイズの分子の RI-MP2 エネルギ
ー計算をルーチン的に実行可能であり、ナノマテリアル設計などへの応用が期待される。
&!#)"

(b)

*+,-./01*+"2345"*."/6467707"*+,-./01*+"8995"*."/6467707"'%("!"##$%"&%!)*"+*!,+*!-$&":76;"*+,"2345"*."/6467707"./012"'!($":76;"*+,"8995"*."/6467707",<067"!"##$%"&

(a)

$&#"

#$%"

$!%&"

#!$%"

!"
!"

#$%"
$&#"
'%("
'%()#*+,-+./0+1,*#&

&!#)"

!"
!"

#!$%"
$!%&"
'!($"
'%()#*+,-+./0+1,*#&

%!)*"図 2. MPI/Open-MP ハイブリッド並列 RI-MP2 アルゴリズムの高速化率(a) タキソール RI-MP2/cc-pVDZ (1123AOs)(b) C60@C60H28 RI-MP2/cc-pVTZ (3992AOs)[参考文献] [1] M. Katouda and S. Nagase, Int. J. Quant. Chem. 109, 2121 (2009).
[2] A. P. Rahalkar, M. Katouda, S. R. Gadre, and S. Nagase, J. Comput. Chem. 31, 2405 (2010).