4E01
E01
GPUを用いた分割統治型電子状態計算の高速化
1（早大先進理工 ，早大理工研2，JST-CREST3, 京大ESICB4）

吉川 武司1，中井 浩巳1-4

Acceleration of divide-and-conquer calculation on GPU(Advanced Science and Engineering, Waseda Univ.1, RISE, Waseda Univ.2, CREST, JST Agency3,
ESICB Kyoto Univ.4)
Takeshi Yoshikawa1, Hiromi Nakai1-4【緒言】
緒言】 近年 GPGPU の量子化学への適用は広範囲に行われている。Hartree-Fock (HF)から高精度な Cupled Cluster (CC)計算にまで幅広く対応しており、大幅な高速化に成功してきた。
その中でも、GPGPU 専用の量子化学プログラムパッケージ『TeraChem』[1]は全てのコードを
CUDA で書かれており、タンパク質等の大規模分子に対しても、HF 計算を数十倍から数百倍の高速化が可能である。しかし、大規模分子に対する高精度な電子相関理論の適用例は少ない。その大きなボトルネックとなっているのが GPU 上のメモリ量である。CPU のメモリに対して GPU 内部のメモリは最大で 5～6 GB と比較的低量である。そのため、大規模分子に対して適用するためには、使用するメモリ量を少なくする必要がある。その解決策として当研究室で開発を行ってきた分割型計算理論の一つである分割統治(DC)法[2,3]の利用を提案する。DC
法等の分割型計算理論の特徴として、必要なメモリサイズを分割することで大幅なメモリ削減ができ、大規模分子への適用が可能となる。本研究では、GPGPU を用いた大規模高精度電子状態理論を可能とするために、DC 法のコードを GPU 化し、高速化を行ったので報告する。
【RI 近似を近似を用いた DC-HF 法の高速化】
高速化】 DC-HF 法は全系をいくつかの部分系に分けて対角化計算を行うことで計算コストを削減する方法である。部分系の周りのバッファ領域を含めた局在化領域で部分系αの軌道を構築することで、周囲の効果を取り込むことが可能となる。
しかし、SCF サイクルにおいて、DC 法は全系の Fock 行列を生成する必要があり、その生成時間が DC-HF 法の大きなボトルネックとなっている。ただし、対角化とは異なり、Fock 行
列の生成は並列化効率が高く、GPU を用いた高速化が可能である。
さらに、
resolution-of-identity[4](RI)近似 を適用することにより Fock 行列生成のコストを下げることも可能である。
DC-RI-HF法ではクーロン項J、交換項KをRI近似より(1), (2)式を用いて求めることができる。
αJ µνDC-RI-HF = ∑ ∑ ∑ ( µν P )( P Q ) ( Q λσ ) Dλσ−1αK µνDC-RI-HF =(1)

λσ ∈S (α ) PQ

−1 2
−1 2
1
α
µσ P )( P R ) ( R Q ) ( Q µσ ) Dλσ(∑∑∑
2 α λσ ∈S (α ) PQR(2)

ここで、µ, ν, λ, σは原子基底を、P, Q は補助基底関数を表す。ここで、J、K を求める部分をGPU で並列に計算を行った。
【 数値検証:
数値検証 DC-RI-HF 法 】 β-ストランド型オリゴグリシン(gly)10 に対して、基底関数は
cc-pVTZ、補助基底関数は cc-pVTZ を用いて DC-RI-HF 計算を行った。DC 法による計算では、部分系はグリシンモノマーからなるユニットとし、バッファは左右 2 ユニットとした。表 1に、direct-SCF 計算と RI-SCF 計算における Fock 行列の構築と対角化の時間を示す。ただし、各項目の高速化率を括弧書きで示す。GPU を用いることによって Fock 行列の構築時間を 4.7

倍程度減少することが可能となる。さらに、RI 近似と GPU を融合させることによって、Fock
行列生成時間を 10.6 倍高速化することに成功した。一方で、並列化が困難である対角化の計
算時間は DC 法を用いることによって高速化することが可能である。今回は Fock 行列生成の部分のみ GPU 化したため対角化の計算時間は GPU を用いた場合でも変化はないが、全てをGPU 化することによって対角化の計算時間も減少することが可能となる。
Table 1. The wall time [min] of construction of Fock matrix anddiagonalization of SCF equation by direct-SCF and RI-SCF on CPU and GPU.
Construction
Diagonalization
Total time
Direct-HF(CPU)
124.2
3.8
128.5
Direct-HF(GPU)
26.7
(4.7)
3.8
(1.0)
30.7
(4.2)
RI-HF(GPU)
11.7(10.6) 3.8
(1.0)
15.9
(8.1)
DC-RI-HF(GPU)
11.7(10.6) 2.5
(1.5)
14.7
(8.8)
CPU: intel Xeon 3.47GHz (1core) GPU: Tesla K20m
【GPU を用いた DC-電子相関法
電子相関法の電子相関法の高速化】
高速化】 DC 法に基づく 2 次 Møller-Plesset (MP2)法や CC法等においては、各部分系の軌道を用いることで電子相関計算を行う。部分系の軌道はバッファ領域にも広がっているため、当研究室で開発されたエネルギー密度解析(EDA)[5]を用いて中央領域だけの相関エネルギーを求めることで、2 重カウントを防いでいる。これを MP2、
CCSD 計算に適用すると、以下のようにして相関エネルギーを見積ることができる。
occ

vir

i, j

a ,b µ∈S (α )DC
∆ Ecorr.
= ∑∑∑ ∑ Cµαi∗ µ j α aα bα  2tɶijα,ab − tɶijα,ba α(3)aα bα iα j α

MP2 の場合:

tɶijα,ab = −CCSD の場合:

tɶijα,ab = tiα,at αj ,b − tiα,bt αj ,a + tijα,abε aα + ε bα − ε iα − ε αj

ここで、ボトルネックとなっている MO 積分、T1、T2amplitude を求める部分を GPU で並列に処理を行うことによって高速化を目指す。

Wall time [hour]【数値検証:
数値検証 DC-CC 法】 オリゴエン C2nH2n+2 に対O(n6.3)して、基底関数は 6-311G**を用いて DC-CCSD 計算を行った。DC 法による計算では、部分系は炭O(n1.4)素二個からなるユニット(C2H4 or C2H5)とし、バッ
15.6x
ファは左右 2 ユニットとした。図 2 に、通常法とDC 法を用いた場合の CCSD 計算に必要な時間とO(n1.3)分子サイズとの関係を示す。ただし、CPU は intel
Xeon 3.47GHz (1core)を、GPU は Tesla K20m を用
Number of units n
いて計算を行っている。DC 法を用いることによFig. 1. System-size dependence of the wallってほぼ線形の計算時間を達成することに成功 times for the DC and conventional CCSD
しているが、CPU を用いた場合のプレファクター calculations of C2nH2n+2 on CPU and GPU.
は依然として大きい。しかし、GPU を用いることによって計算時間の大幅な削減に成功し、15.6 倍の高速化に成功した。
[1] I. S. Ufimtsev and T. J. Martinez, J. Chem. Theor. Comput. 4, 222 (2008).
[2] W. Yang, Phys. Rev. Lett. 66, 1438 (1991).
[3] M. Kobayashi and H. Nakai, in Linear-Scaling Techniques in Computational Chemistry and Physics: Methods andApplications (2011, Springer), pp. 97-127.
[4] J. L. Whitten, J. Chem. Phys. 58 (1973) 4496. [5] H. Nakai, Chem. Phys. Lett. 73 (2002) 363.

4E02
大規模・高精度相対論的量子化学理論の開発（４）：内殻固定ポテンシャル法の開発（早大先進理工 1、早大理工研 2、JST-CREST3、京大 ESICB4）○清野淳司 1、樽見望都 1、中井浩巳 1,2,3,4

Large-scale and highly accurate relativistic quantum-chemical theory (4):Development of frozen core potential method(Advanced Science and Engineering, Waseda Univ.1, RISE, Waseda Univ.2, JST-CREST3, ESICB, Kyoto Univ.4)○Junji Seino1, Moto Tarumi1, and Hiromi Nakai1,2,3,4【緒言】我々はこれまで 4 成分 Dirac 法と同等の精度で電子状態の記述が可能な 2 成分相対論、多電子系のための無限次 Douglas-Kroll-Hess（IODKH）変換法[1,2]を基盤に、局所ユニタリー変換（LUT）法の開発[3,4]および分割統治（DC）法への拡張[5]を行うことで、大規模分子系のための相対論的量子化学理論の開発を行ってきた。この手法は十分な相対論効果および電子相関の双方を含み、かつ計算全体の線形スケーリングを達成する方法である。しかし、重元素を多数含んだ
大規模分子系を計算する場合、あらわに取り扱う電子数の多さは計算の大きなボトルネックとなる。従来、重原子化合物を扱う多くの研究では、相対論的な内殻ポテンシャルを用いた有効内殻
ポテンシャル（ECP）法やモデルポテンシャル（MP）法が使われてきた。これらの手法は分子構造や価電子が支配的な分子物性計算において多大な成功を収めてきたが、重元素が直接関与する
化学反応や電磁気的物性などにおいては多くの議論の余地がある。そこで本研究では、ECP や
MP 法のようなポテンシャル法と相対論的な全電子計算をスムーズに結合する手法、内殻固定ポ
テンシャル（FCP）法を提案する。本手法はポテンシャル法の限界を精査できる可能性を有すると同時に、我々が開発してきた LUT-IODKH 法の大幅な高速化が見込める方法でもある。
【理論】FCP 法は MP 法で用いられる Huzinaga-Cantu 方程式[6]を基盤に構築される。FCP の価電子ハミルトニアンは次のように書ける。
H FCP

v

nv


i

h 2 (i) nv

g 2 (i, j ) 

i j

Q Q

 r




12 

nv

 V (i)
i

λ
FCP

(1)



ここで、h 2 、g 
2 はそれぞれ 1 電子・2 電子ハミルトニアンを表し、Q は有効核電荷である。また、
λは原子で決定された内殻ポテンシャルであり、次のように定義される。
VFCP

J  K  P

λ
VFCP
2

λ
c

c

λ
c

λ

(2)

c

ここで、 J cλ は内殻のクーロン演算子、 K cλ は内殻の交換演算子であり、 P λ は価電子が内殻軌道に落ちることを防ぐための射影/シフト演算子である。通常 MP 法では、h 2 内の一部の核―電子相互作用、および J cλ 、 K cλ が原子計算で決定されたパラメータを伴った関数で表現されるが、本手法では原子計算で得られた内殻軌道をそのまま用いる。また、h 2 や g 
2 は IODKH 法などによる電子状態ハミルトニアンを直接用いることで、高精度な相対論効果を取り入れることが可能である。
特に 1 電子・2 電子共に IODKH 法を適用することで、電子で閉じた（陽電子の情報を一切含まない）最も高精度な相対論的価電子ハミルトニアンを与える。
【結果と考察】FCP 法による数値検証をするため、貨幣金属 2 量体（Cu2、Ag2、Au2）の計算を行った。ハミルトニアンの 1 電子項として非相対論（NR）、1 次から 3 次の DKH 法（DKH1、DKH2、
DKH3）、IODKH 法とその LUT 法（IODKH、LUT-IODKH）を、2 電子項として非相対論的クーロ

ン相互作用（NR）と LUT 法を適用した IODKH 法による相互作用を用いた。また、基底関数として NR/NR のみ Sapporo-(DKH3)-DZP-2013[7]を、その他は DKH3-Gen-TK/NOSec-V-TZP [8]を用いた。Table 1 に全電子計算（AE）による全 SCF エネルギーと FCP の AE からの差を示す。この結果、FCP の差、つまり固定された内殻軌道の分子形成によるエネルギー誤差は、大部分で 0.00001
hartree と非常に小さいことがわかる。このように分子を形成に伴う、内殻軌道の変化に対するエ
ネルギー寄与が非常に小さいため、ECP や MP 法は分子構造などで非常に有効であると推察される。また、それぞれの 2 成分相対論における価電子エネルギーの精度を検証するため、Table 2 に
LUT-IODKH/IODKH による価電子エネルギーとその他の手法の差を示す。NR/NR の差は価電子エネルギーの全相対論効果に対応し、(LUT-)IODKH/NR の差は 2 電子の相対論効果に対応する。この結果、2 電子の相対論効果の全相対論効果に占める割合は非常に小さいが、効果自体は数十
millihartree と無視できないことがわかる。これにより価電子においても 2 電子の相互作用が分子物性に影響を及ぼす可能性があることが示唆される。
Table 1: Total SCF energies of AE and its deviations of FCP from AE (in hartree).
1e Hamil.

LUT-IODKH/
NR/NR

DKH1/NR

DKH2/NRIODKH/NR/2e Hamil.
AE

LUT-IODKH

DKH3/NRNR/IODKH-3277.863087 -3312.473032 -3306.252807 -3306.396210 -3306.389474 -3306.389466 -3306.879946

Cu2
FCP*
AE

0.000007-0.000002

0.000008-0.000003

0.000007

0.000007-0.000003-10395.209983 -10680.451017 -10623.781542 -10625.979712 -10625.824207 -10625.824210 -10628.298985

Ag2
FCP*
AE

0.000009

0.000071

0.000012-0.000001-0.000005-0.000004

0.000008-33089.179424 -38645.529625 -37985.652959 -38027.221489 -38022.392886 -38022.392874 -38037.816243

Au2
FCP*

0.029665

0.000069-0.000072-0.000016-0.000023-0.000023-0.000019

Table 2: Valence SCF energies of LUT-IODKH/IODKH and its deviations of the others from LUT-IODKH/IODKH in FCP (in hartree).
1e Hamil.

LUT-IODKH
NR/NR

DKH1/NR

DKH2/NR

LUT-IODKH

DKH3/NRIODKH/NR/2e Hamil.

/NR/IODKHCu2

1.205092-0.027350

0.025983

0.024733

0.024828

0.024833-392.879709

Ag22.328446-0.095973

0.036324

0.031100

0.031272

0.031273-293.065054

Au2

10.189726-0.550251

0.105630

0.064195

0.066943

0.066944-272.328174

また、FCP 法の価電子エネルギー計算において、内殻関数を削除した基底関数を用いたときの計算時間を Table 3 に示す。削除した関数は Cu、Ag、Au でそれぞれ[2s1p]、[3s2p]、[4s3p2d1f]である。

この結果、FCP 法は全電子計算に比べて最大 51.1 倍の高速化が達成された。さらに大きな分子系ではより有効であると予想されるため、大規模重原子化合物
計算において非常に有用な手法であることが示された。

Table 3. CPU time (in second) in FCP with
reduced basis sets and AE calculations.

【参考文献】[1] M. Barysz and A. J. Sadlej, J. Chem. Phys. 116, 2696
Mole.
(2002). [2] J. Seino and M. Hada, Chem. Phys. Lett. 461, 327 (2008).
[3] J. Seino and H. Nakai, J. Chem. Phys. 136, 244108 (2012). [4] J.
Cu2Seino and H. Nakai, J. Chem. Phys. 137, 144101 (2012). [5] J. Seinoand H. Nakai, J. Chem. Phys. 139, 034109 (2013). [6] S. Huzinaga andAg2
A. A. Cantu, J. Chem. Phys. 55, 5543 (1971). [7] T. Noro, M. Sekiya,
Au2
and T. Koga, Theor. Chem. Acc. 132 1363 (2013). [8] Y. Osanai, T. Noro,
E. Miyoshi, M. Sekiya, and T. Koga, J. Chem. Phys. 120, 6408 (2004).
*t(AE) / t(FCP)

t

FCP

AE

Speed-up*

t

6.0

(2.2)

13.2

19.1

(4.5)

85.5

33.2

(51.1)

1698.14E03
2 成分相対論法を用いた同位体分別反応における核の体積効果の理論的研究（首都大院理工*、早大先進理工**、JST-CREST***）○根本佳介*,***、阿部穣里*,***、清野淳司**,***、波田雅彦*,***An ab initio study of nuclear volume effects using 2-compornent
relativistic method for isotope fractionation reaction（Tokyo Metropolitan Univ. *, Waseda Univ. **, JST-CREST***）○Keisuke Nemoto*,***, Minori Abe*,***, Junji Seino**,***, Masahiko Hada*,***【諸言】鉛やウランなどの重原子の同位体組成の分析は、年代測定など地球化学・分析化学の分野で広く用いられている。同位体組成は化学平衡過程をへて変化しており、この様な同位体分別は自然界でも広く確認されている。一般的な同位体分別反応 A’Y + AX = AY + A’X (A’ : 軽い同位体、A : 重い同位体) において、平衡定数はほぼ 1 に近く、平衡定数から 1 を引いた値を同位体分別係数と呼ぶ。平衡反応における同位体効果は、同位体原子核の質量の違いに起因する振動準位の違いを考慮する理論として Bigeleisen らによって説明され 1、約半世紀の間、同位体科学の分野で支持されてきた。しかし、重原子の同位体分別においては、分子振動の効果だけでは説明することができなかった。1996 年に藤井らの実験 2 及び Bigeleisen の理論 3 により、重原子の同位体分別は主に原子核の体積差により生じる電子状態の差による影響であることが明らかとなった。
核体積効果は重原子では顕著となり、また核近傍の電子状態が影響するため相対論効果が極めて重要となる。そこで Schauble4、阿部ら 5,6 は同位体分別係数における核体積項を見積もる量子化学的手法を、4 成分相対論法に基づいて提案した。4 成分相対論法は相対論的に大変正確であるが計算コストが高く、数原子分子への適応に限られる。より現実的な分子系への適応に向け、精度を落とさずにより大きな分子が計算できる理論の選定が望まれる。本研究では、4 成分相対論法の近似的な相対論法として有力な、Douglas-Kroll (DK) 法に基づく 2 成分相対論法を用いて、核の体積効果が顕著であるウラン分子で検証を行った。スピン-軌道相互作用や高次の相対論効果の影響について、4 成分相対論法との比較をHartree-Fock 法を用いて行った。
【理論・計算】核の体積効果の理論計算における相対論効果の影響を調べるために、ウラン化合物を用いた同位体分別として、4 価ウラン (UCl4) － 6 価ウラン
(UO2Cl3-, UO2Cl42-) 系における 238U、235U の同位体交換反応 235U (Ⅳ) + 238U (Ⅵ) = 238U (Ⅳ) + 235U (Ⅵ)
を考えた。ここで同位体分別係数における核の体積

Fig. 1 ウラン４価－6 価の反応物モデル

項 lnKnv は式 (1) で表される。
-1

lnKnv =(kT) {[E (238

U(Ⅵ)) -E (235

U(Ⅵ))] - [E (238

U(Ⅳ)) -E (235

U(Ⅳ))]}(1)

U 同位体の核半径は文献値 (U238 : 5.8514 [fm], U235 : 5.8236 [fm]) を用い、Gauss 型有限核モデルに使用した。
計算方法について Table 1 にまとめた。
計算プログラムは DIRAC08 及び DIRAC10を使用し、基底関数は U : (25s21p16d12)、Cl : (16s11p2d)、O : (14s9p2d)を用いた。4 成分相対論法として Dirac-Coulomb 法 (DC)、2 成分相対論法として DK 法を用いた。より高次の DK 法を用いることで計算精度は高くなり、中でも Infinite-Order DK (IODK) 法は 1 電子ハミルトニ
アンに関して 4 成分 Dirac ハミルトニアンと等価となる。そこで、相対論効果における高次の DK
変換項の影響を非相対論法、2nd-Order DK 法 (DK2) 、IODK 法にて比較を行い、スピン-軌道
相互作用の影響についても spin-free (SF)、1 次 (SO1)、2 次 (SO2)、無限次 (IOSO)ついて確認した。また 2 電子演算子のスピン-軌道相互作用を近似的に取り込む方法として、原子平均場に基づく spin-same-orbit 項と spin-other-orbit 項を考慮した Mean-Field Spin-Orbit (MFSO2) 法の影響も評価した。
Table 1 ハミルトニアンの分類
非相対論
Method
(NR)
SF
SF
0
2次
SO
0
0
2電子SO
0
0

DK2
SO1
2次
1次
0

IODK
4成分相対論(DC)SF
SO1
IOSO IOSO+MFSO2
無限次 無限次 無限次
無限次
無限次
0
1次
無限次
無限次
無限次
0
0
0
MF
無限次

SO2
2次
2次
0

【結果及び考察】ウラン分子系における体積効果 lnKnv 値の計算結果を Table 2 に示した。スピ
ン非依存の結果として、NR、DK2/SF、IODK/SF、DC の lnKnv 値を比較すると、高次の DK 変
換項の寄与は無視できず、非相対論法や低次の DK 法では正しく計算できないことがわかる。しかし IODK/SF においても DC に対して約 1‰程度の差がみられた。スピン-軌道相互作用の影響について IODK/SF, SO1, IOSO での lnKnv 値を比較すると、高次項についても効果が大きく、SO1
では IOSO と比べ lnKnv 値が高く見積もられることがわかる。また IODK/IOSO に 2 電子演算子のスピン-軌道相互作用として、MFSO2 法を加えることで、DC とほぼ同精度となった。そのた
め 2 電子演算子のスピン-軌道相互作用の影響も lnKnv 値を計算する上で重要であることがわかった。計算時間は UCl4 分子のエネルギー計算の一点計算において、4 成分相対論法: 84 hour、
IODK/IOSO + MFSO2: 1 hour 17 min であった。IODK 法にスピン-軌道相互作用を適切に考慮した理論 (MFSO2) を用いることで 4 成分相対論法とほぼ同精度になり、計算コストを大幅に削減できることを明らかにした。今後はより現実的な分子系への応用として、50 原子程度の中規模分子へ核体積効果の理論計算を用いた解析に適応できると考えられる。
Table 2 各計算レベルにおける核の体積項(lnK nv [‰])計算結果 (T = 308K)
MethodSO[UCl4] - [UO2Cl4]2[UCl4] - [UO2Cl3]1-

0.03

SF
1.30

DK2
SO1
2.12

SO2
1.88

SF
1.57

SO1
2.50

IODK
IOSO IOSO+MFSO2
2.31
2.51

2.53

-0.28

1.24

2.09

1.88

1.61

2.48

2.30

2.53

NR

2.47

DC

Expt.[7](U[Ⅳ]-U[Ⅵ])2.24【参考文献】[1] J. Bigeleisen et. al., J. Chem. Phys., 15, 261 (1947). [2] J. Bigeleisen, J. Am.

Chem. Soc., 118, 3676 (1996). [3] M. Nomura et. al., J. Am. Chem. Soc., 118, 9127 (1996). [4] E.
A. Schauble, Geochim. Cosmochim. Acta, 71, 2170 (2007). [5] M. Abe et. al., J. Chem. Phys.,
128, 144309 (2008). [6] M. Abe et. al., J. Chem. Phys., 133, 044309 (2010). [7] Y. Fujii et. al., J.

Nucl. Sci. Technol., 43, 400 (2006).

相対論的分子軌道法における種々の２成分法の近似精度について（九大院理） ○井上 頌基，鈴木 聡，渡邉 祥弘，中野 晴之

Accuracy of the two-component methods in the relativistic molecular orbital theory(Kyushu Univ.)○Nobuki Inoue, Satoshi Suzuki, Yoshihiro Watanabe, and Haruyuki Nakano
ion@ccl.scc.kyushu-u.ac.jp【序】原子分子における相対論効果を 4 成分法に比べて低コストに取り扱える 2 成分法は，近年までに数多くの方法が提案されている．これらの 2 成分法がそれぞれどの程度の精度を有するか検討することは計算方法の選択のために有意義であると考えられる．そこで本研究では代表的な 2 成分法（Breit-Pauli 近似; BPA [1], ZORA[2], IORA[3], RESC[4], DK[5], IOTC[6]）のふるまいと近似精度を 1 電子演算子と 2 電子演算子に分けて議論する．
【方法】 1 電子系および 2 電子系に対し，BPA,ZORA,IORA,RESC,DK,IOTC および非相対論（NR），4 成分法（4-comp.）それぞれに関して，プログラムの実装を行い，算出した物理量を比較した．ここで，2 電子反発エネルギーはクーロン反発項（g (i ,j) = r ）として計算したが，特に精度の高い DK と IOTC については，それぞれ DK 変換，IOTC 変換を施した 2 電子反発演算子（g (i ,j) , g (i ,j)）も実装し，計算を行った．g (i ,j)は[7]の論文において，
既に導出がなされている．しかし，g (i ,j)については g (i ,j)以外の高次の場合において，平均場ではなく演算子としての具体的な表式は知られていない．そこで，g (i ,j)，g (i ,j)に関して，その具体的な表式を新たに導出した．

E : atomic unit

【結果】 ①１電子系 核電荷 Z=96 のクーロン場に対し，単一の s 型ガウス型波動関数をもつ
１電子のエネルギーE を計算した．E とガウス関数の指数との関係を図.1 に示す．縦軸は Eを，負の値でも扱えるようにした対数であり，横軸は指数の対数である．このグラフから，
2 成分法によって見積もられるエネル
NR
ZORA
ギーについて以下の知見が得られる．
(i) 第一は，が大きいときの漸近的な挙動である． は運動量の大きさの２4-comp.
乗（p2）に比例するため，と E の絶対値が十分に大きな領域では，グラフの傾きから p2 に対し（運動）エネルギーがおおよそ何次に比例するかを読み取れる．
(ii) 第 二 は ， BPA 以 外 は す べ て =10000 付近に最小値をもつことから，変分的に軌道を決定するための必要条
件は満たしていることである．
BPA(iii) 第三は，近似精度である．プロッ
トしたの全区間において DK は非常に精度良く 4 成分法を近似しており，
log10 また，IOTC は完全に 4 成分法のエネル
ギーを再現しているため，
図.1 1 電子１中心系における s 型ガウス関数のグラフ上ではこれらは 4 成分法の線とエネルギーに対する指数依存性
重なっている．

E : atomic unit

②２電子系における 電子ハミルト
ニアンの精度 核電荷 Z=96 のクーロン場において，単一の s 型ガウス型波動関数をもつ 2 電子系の全エネルギ
ーE を計算した．1 電子ハミルトニア
ンはおのおのの 2 成分法で，2 電子反発項は非相対論で評価している．E
とガウス関数の指数  の関係を図.2に示す．図.1 と図.2 を比較すると非常によく似ていることから，全エネ
ルギーの精度に大きく寄与するのは１電子ハミルトニアンの近似精度であることが示唆される．

log |(E-E)/E|

E : atomic unit

③２電子系における２電子反発演算
log10 子まで含めた近似精度
DK2 と DK3 について 2 電子反発演
図.2 2 電子１中心系における s 型ガウス関数の算子を導出した． DK1，DK2，DK3，
エネルギーの指数依存性
IOTC について②と同様の系に対して全エネルギーE を計算した．ここでは１電子ハミルトニアンを DK1と IOTC で評価した場合の E の 4 成
分法との相対値とガウス関数の指数との関係を図.3 に示す． “/”で区切って１電子ハミルトニアンの計算
方法と２電子反発項の計算方法を示している．DK1 の１電子ハミルトニ
アンの誤差が２電子項の誤差に比べて大きいため，DK1/NR と DK1/DK1の誤差はほぼ一致している．一方，
IOTC の１電子ハミルトニアンは十分に誤差が小さく，２電子演算子のlog10 log10 
精度が高いほど，誤差が小さくなっていることが分かる．最も精度が高
図.3 2 電子１中心系における s 型ガウス関数のエネルギーの指数依存性
い IOTC/IOTC では 1 電子系のように4 成分法と完全には一致しないものの，おおよそ６桁程度の精度が期待できる．今回導出した式に基づく IOTC/DK2 と IOTC/DK3は，IOTC/DK1 と IOTC/IOTC の中間的な値をとっており，近似の次数から予測される通りの結果である．
結果のより詳細な解析，および，その他の結果については当日報告する．
[1] Bethe, H. A.; Salpeter, E. E. Quantum
mechanics of one- and two-electron atoms; Springer: Berlin, Heidelberg, New York, 1957.
[2] E.van Lenthe, E.J.Baerends, and J.G.Snijders, J. Chem. Phys. 101, 9783 (1994)[3] K. G. Dyall and E. van Lenthe, J. Chem. Phys.111, 1366 (1999)[4] T. Nakajima and K. Hirao. Chem. Phys. Lett.,302,383-391 (1999)[5] T. Nakajima and K. Hirao. Chem. Phys. Lett. 329 511-516 (2000)[6] M. Barysz and A. J. Sadlej, J. Chem. Phys. 116, 2696 (2002)[7] J.Seino and M.Hada. Chem. Phys. Lett., 461, 327-331 (2008)

4E05
軌道エネルギーにもとづく化学反応経路解析法（山梨大燃研） ○常田 貴夫
Chemical reaction path analysis based on orbital energies(Univ. Yamanashi)○ Takao Tsuneda【序】 化学反応解析は，反応過程の反応障壁にもとづいて行なわれ，分子の電子状態をあらわに考慮した定量的な議論は通常なされない。特に気相反応の反応性は反応段階での反応障壁エネルギーで議論されるのが主であり，まれに状態密度や反応断面積が考慮される。Frontier 軌道
論は気相反応の反応性を分子の電子状態間の相互作用として議論するまれな例である[1]。この方法では，分子軌道と対応する軌道エネルギーを使ったダイアグラムで反応性を議論する。しかし，分子軌道の位相のみによる定性的な議論が主であり，軌道エネルギーは参考程度の利用にとどまる。これは，軌道エネルギーを定量的に再現できる理論がなかったことに起因する。
現在，軌道エネルギーを定量的に再現できる理論が存在する。長距離補正（LC）密度汎関数
法（DFT）[2]である。LC は，DFT 交換汎関数が考慮していない長距離交換を，HF 交換積分の長距離部分と組み合わせて補完する補正である。LC はその簡便さにもかかわらず，DFT 計算の問題とされてきた van der Waals 結合や電子移動励起，超分極率など非線形光学応答物性が再現できない問題を包括的に解決する[3]。最近，この LC 汎関数を利用した Kohn-Sham(KS)法は，
B3LYP など従来の汎関数を使うと大きく過小評価される価電子軌道エネルギーを 0.1eV 程度の
平均絶対誤差で定量的に再現することがわかった[4]。一般的に，占有・非占有軌道エネルギーはイオン化ポテンシャルと電子親和力の逆符号にそれぞれ一致する[3]。この高精度な軌道エネ
ルギーを利用すれば，軌道エネルギーにもとづく定量的な反応解析が可能になるはずである。
軌道エネルギーと化学反応との関係性について，最近興味深い結果が得られた。Diels-Alder
反応の固有反応座標（IRC）に沿って LC-KS 計算による HOMO-LUMO ギャップをプロットした結果，反応初期段階でほぼ一定で，遷移状態付近で急増することがわかった[5]。IRC 上の軌道エ
ネルギーギャップの傾きは，構造変化と電子移動にともなう変化である。電子数が増減する軌道の軌道エネルギー変化はゼロである[3]ため，反応初期段階でのほぼゼロの軌道エネルギーギャ
ップの傾きは，反応が構造変化の影響なしに電子移動で進行する高速反応であることを意味する。
これを利用すれば，対象の化学反応経路が高速反応かどうかを解析することができる。
本研究では，IRC 上の反応に関与する軌道の軌道エネルギーギャップの傾きを利用し，軌道エネルギーにもとづく化学反応経路の反応性評価法を提案する[6]。

【方法】反応性評価は次の手順で行なう。
(1) LC-DFT 計算により，IRC を求めて価電子エネルギーをプロットする。
(2) 軌道と軌道エネルギーの変化をもとに，価電子軌道で最も安定化・不安定化する占有・非占有軌道を「対象軌道」として選び出す。
(3) 対象軌道の軌道エネルギーギャップを
IRC 上にプロットして図１の規格化反応
ダイアグラムを作成する。気相反応では反応物と生成物の錯合体を新たに反応物と生成物と定義する。錯合体を形成しない場合は，電子移動が起こる
直前の状態を反応物と生成物とする。
(4) 反応初期段階での軌道エネルギーギ
ャップの傾きを計算して反応性指標と図１．規格化された反応ダイアグラムの例。

し、その絶対値が 0.25 以下なら電子移動のみにより反応が進行していると解釈する。
以上の方法で求めた反応性指標を利用し，化学反応を解析した。

【結果】43 種類の基本的な化学反応について反応ダイアグラムを解析した結果，6 割以上について順反応の初期段階が電子移動のみで進行することが示された。一方逆反応では 4 割弱と少なく，
軌道エネルギーギャップが小さくなる反応性が低い反応が多いことが分かった。さらに，反応障壁の高さが順反応と逆反応で逆転している反応がいくつか存在するため，これらの反応の反応性指標を比較した結果，1 反応を除いて逆反応が電子移動で進行していないことが示された。その唯一の例外である CH4 + Cl→CH3 + HCl 反応については，調査の結果，IRC とはまったく異なる反応経路を通ることが衝突実験で示されていることがわかった[7]。反応ダイアグラムから予想される反応性と矛盾する反応も存在する。これらの反応はおおむね SN2 反応，対称反応，CH3 反応，そして開環反応の 4 種類に分類できる。まず SN2 反応については，交差分子ビームによる衝突実験により，IRC とはまったく違う反応経路を通ることが実証されていることがわかった[8]。すなわち，
SN2 反応の Cl- + CH3I → CH3Cl + I-反応の時間分解解析の結果，この反応は迂回経路，つまりCH3 基が I-のまわりを回転しながら Cl-に移動する反応経路をとる。他の SN2 反応も迂回経路をとるものが存在すると予想できる。対称反応は反応物と生成物が同じ反応のことを指し，主な反応は水素原子の交代か移動のどちらかである。対称反応はほとんど実験的に解析されていないが，水素原子交代による対称反応（たとえば，H + HCl→HCl + H）は SN2 反応によく似ており，類似した反応経路で進行する可能性がある。また，水素原子移動による対称反応（たとえば，pentadieneの二重結合交代）については，きわめて大きなトンネル効果が予想され，本研究のようにトンネル効果を考慮せずに反応を議論するのは難しい。CH3 反応については，CH4 + Cl→CH3 + HCl 反応
以外の反応は順反応の反応障壁は逆反応より圧倒的に低い。詳細は当日の発表で紹介する。

【参照文献】[1] K. Fukui et al., J. Chem. Phys. 20, 722, 1952. [2] H. Iikura et al., J. Chem. Phys.
115, 3540, 2001. [3] 常田貴夫, 「密度汎関数法の基礎」(講談社, 2012). [4] T. Tsuneda et al., J.

Chem. Phys. 133, 174101, 2010. [5] R. K. Singh and T. Tsuneda, J. Comput. Chem. 34, 379, 2013.
[6] T. Tsuneda and R. K. Singh, submitted. [7] W. R. Simpson et al., J. Chem. Phys. 103, 7313, 1995.
[8] J. Mikosch et al., Science, 319, 183, 2008.

4E06
From Atom to Molecule 法と increasing exchange 理論による有機化合物のシュレーディンガー解の計算(量子化学研究協会・JST-CREST)

中辻

博、中嶋

浩之

From-atom-to-molecule idea combined with increasing-exchange theory:
Schrödinger solutions for organic molecules(QCRI, JST-CREST) Hiroshi Nakatsuji, Hiroyuki Nakashima
化学の中で最も成功している概念は、「化学構造式」であろう。このことは、分子は原子から３次元的に構成されているという概念が普遍的に正しいことを示している。我々はこの概念を分子の波動関数の構成に利用し、From Atom to Molecule (FATM) 法を提案した。この理論を有機化学に応用するにあたり、Pauli 原理に基づく反対称化を、原子内電子の反対称化と原子間のそれとに分けて行う increasingexchange (iExg) theory によって表現し、大きな簡単化と加速を実現した。この FATM-iExg 理論は化学的直観に富んだ自然な order-N 理論であり、有機化合物のシュレーディンガー解の計算に適している。

シュレーディンガー方程式の正確な解法[1,2]に従えば、Ne 電子系の正確な波動関数 は、
対 象 と す る 分 子 の ハ ミ ル ト ニ ア ン と ス レ ー タ ー 型 VB 近 似 波 動 関 数  0 か ら 、 Free
Complement (FC)法によって自動的に生成されψ   cI  I(1)Iと書かれる。この方法を「直接法」と呼ぶ。 I は完員関数と呼ばれ、次式で与えられる。
Ne
Ne
NA


 kr ,iA k mij kkk
 I 1, 2..N     exp( i riAi ) xiAxi,iAi yiAyi ,iAi ziAzi,iAi riAi i   riA iA   rij i 
A (  Ai )j ( i )(2)

完員関数は Hartree 積の形で書かれ、i,j,e は電子を、A は原子を、Ai は電子 i が属する原子を示す。電子の原子への所属の仕方は自由でありイオン項も現れる。冪にある k や m は完員関数の次数に関係する整数(0,1,2..)である。(2)式の完員関数には、原子部分の積のみからなるものと、原子間にまたがる関数からなるものがある。
この関数のハミルトニアン積分や重なり積分の計算は一般の原子・分子では不可能である。
そのため変分法を使うことはできず、分子座標上の任意のサンプリング点でシュレーディンガー方程式が成立しているという必要条件、Local Schrödinger equation (LSE) 法によって計算してきた[3]。その結果はレビュー[1]等をご覧いただきたい。
これに対して、ここで述べる FATM(from atom to molecule)法では、この波動関数をより化学的に理解しやすい化学構造式の原理に基づいて構成する。結果は直接法と等価である。
化学構造式は化学の日常で使われている大変有用な原理であり、分子は原子の３次元的な配列から成り、化学結合もそれによって特徴づけられる事を示している。この考えを基礎に、
例えば２原子分子 AB の正確な波動関数を     cijXY A XY  iX Yj   cIXY A XY inter
XY , I 
I( XY )  ij
(3)で表現する。ここで XY は Ne 電子の核 A,B への所属の自由度から生じる共有結合項 AB、イオン項 A+B-、A-B+など異なる電子配置項を示し、最初の和はその異なる電子配置項についての和を示す。第一項は原子波動関数とその間の交換相互作用からなる項で、 i は原子 X の i
X番目の状態である。第二項は２中心の関数からなるもので、原子間の電子状態を詳述する部分である。一般の分子では３中心・４中心項・・も現れる。Heitler-London の研究からも分かる通り、化学結合の大きな部分は第一項で書かれる。イオン項の存在も重要である。状態 i, j
を多くとればこの第一項だけで分子を記述することも可能である。この第一項の  iX , Yj を原子 X と Y の完員関数 iX , Yj で書き、その係数を分子で reoptimize することで電子状態の再構成を書くこともできる。(1-2)式と比較して各項の対応は自ずと明らかである。 inter
XY , I は(2)式で定義される原子間にまたがる完員関数であり、それほど高次項は要らない。
(3)式の A は反対称化演算子であり、計算的にも時間がかかるので、効率の高い理論が必要である。関数がスレーター行列で書ける時には N 3 / 3 の速さであり妥当であるが、 rij 項が複雑に入ると遅くなる。我々はすでに FC 型の波動関数について効率の良い反対称化アルゴリ
ズムを提唱しており[4]、普通の分子では  N 35 のオーダーである。ここでは FATM 法に適した反対称化法である漸増交換(increasing exchange (iExg))法を使う。この理論は 2010 年 6月に著者の一人により提唱されていたが、最近実用化された。
上の分子 AB の一つの電子配列項 XY を例にとると、系の反対称化演算子 A XY は次のように書かれる。
(1)
(2)
(3)
(K )
A XY  A X AY (1  E XY
 E XY E XY   E XY
)

(4)

(i )

ここで A X は原子 X に属する N X 個の電子の反対称化演算子、 E XY は X・Y 間の i 対の電子交
換演算子、K は N X と NY の小さい方の数である。ただし N X  NY  N e 。XY 間の一対の電子の交換ごとに関係する XY 間の軌道の重なり積分の２乗の確率が掛る。従って(4)式のうち高次の交換のエネルギーへの寄与は事実上ゼロになる。このことから、総数として ( N X  NY )!
個ある反対称化操作も実際は(4)式の最初の数項で十分になってしまう。(4)式を多原子分子に拡張することは容易であり、時間のかかる反対称化操作は、交換する原子間の距離の増加と共に消滅してしまい、計算は容易にオーダーN になる。有機分子の場合、炭素は６電子で、その反対称化は 6!  720 個と小さく、原子間も隣接でも高々３価電子交換どまり、離れると指数関数的にその数は減り、C60 といえども多くは減衰してしまって実現可能な数にとどまる。
我々はこの FATM－iExg 法によって、まずは炭化水素、次に N を含む有機化合物などに、
シュレーディンガー精度の量子化学を確立していきたいと考えている。
[1] H. Nakatsuji, Acc. Chem. Res. 45, 1480 (2012). 中辻 博, 最先端研究:「シュレーディンガー精度のシミュレーションを目指して」,シミュレーション,

32, 39 (2013). [2] H. Nakatsuji, Phys. Rev. Lett. 93,

030403 (2004). [3] H. Nakatsuji, H. Nakashima, Y. Kurokawa, A. Ishikawa, Phys. Rev. Lett. 99, 240402(2007). [4] H. Nakashima, H. Nakatsuji, J. Chem. Phys. 139, 044112 (2013).

4E07
FATM-iExg 法によるシュレーディンガー解の計算: 炭化水素系有機化合物への応用(量子化学研究協会研究所, JST-CREST) 中嶋 浩之, 中辻 博

Solving the Schrödinger equation with the FATM-iExg method: Application to hydrocarbons(Quantum Chemistry Research Institute (QCRI), JST-CREST) Hiroyuki Nakashima, Hiroshi Nakatsuji
1. FATM-iExg 法Free Complement - Local Schrödinger Equation (FC-LSE) 法は、原子・分子の精密なシュレーディ
ンガー解を求める方法として提案され、これまでに様々な系に適用された[1,2]。系を定義するハ
ミルトニアンが、正確な波動関数の関数空間: 完員関数(Complement function)を自動的に生成し、その未知係数は、積分計算を要する変分法や積分フリーな Local Schrödinger Equation (LSE) 法により決定される。また、「化学式」に代表される「Locality」と「Transferability」を利用し、分子の波動関数の効率的な記述を可能にする From Atom to Molecule (FATM) 法と、電子の反対称化法の新しい理論である Increasing Exchange (iExg) 法が提案された。FATM 法では、原子核近傍の記
述に原子の精密波動関数を再構成して使うことで、結合領域の波動関数は低オーダーで済ませられる。iExg 法では、反対称化演算全体を原子内の電子の反対称化と原子間のそれに分け、後者は原子間距離が増すほど指数関数的に減衰するという原理を利用する。そのため、理論上大きな分子系で自然な Order-N 理論を実現することができる。これらの理論の詳細は、先の講演: 4E06 (中辻, 中嶋)にて述べられる。
我々は、まず有機化合物のシュレーディンガー解の計算技術の確立を目標とし、本稿では主に炭化水素系有機化合物に適用し、理論の検証と安定な計算手法の開発を行った。これらの分子は炭素と水素のみから成るため、まず、炭素原子の精密な計算と(水素原子は自明)、C-H 結合、C-C結合の正確で効率的な記述法を確立することが大切で、これによりあらゆる炭化水素系有機化合物の計算が可能となる。そのため、原子の精密な計算、LSE 法での安定な計算法(サンプリング点の取り方と収束法)、結合領域の効率的な記述、を系統的に詳細に検討した。
2. 原子の計算: LSE 法での安定な計算法(サンプリング点と収束法)
積分フリーな LSE 法は原理上どのような原子・分子にも適用することができる。しかし、サン
プリング点依存性や数値的不安定性を軽減する工夫が必要である。特に原子の計算では、サンプ
リング点が互いに近距離であり、サンプリング点の取り方に注意が要る。また、波動関数の記述が不十分な領域のローカルエネルギーは大きく振動する。このようなサンプリング点依存性を軽減するため、H-square error を指標に記 表 1. C 原子(n=2, Mn=368)の計算
繰り返
し

サンプリン
グ点

Energy (a.u.)

ΔE (a.u.)

H-squareerror

点の H-square error を悪くする点が取り

0
1
2
3
4

3000000
2999502
29991122999063
2999049-37.836 34-37.842 73-37.845 72-37.845 17-37.845 23

8.65  10-3
2.26  10-3
-7.20  10-4
-1.75  10-4
-2.34  10-4

4.49  10-1
3.65  10-1
3.57  10-1
3.57  10-1
3.57  10-1除かれた。その結果、n=2 であっても絶Exact

述の悪い点を取り除く収束法が提案された。表 1 に炭素原子(n=2, Mn=368)の計
算結果を示す。収束法を適用したところ、
全 300 万のサンプリング点のうち 951-37.845

対エネルギーとして 0.1 mH の誤差のオーダーまで正しい解を得ることができた。大きな分子の計算のために、コンパクトでかつ精度の良い原子の精密波動関数を得ることができた。
3. 結合領域の記述

表 2. FATM-iExg 法による CH 分子の計算2 で得た原子の精密波動関数を用い、C 原子: n=2, 原子間: n=1

FATM 法に基づき分子 C 原子: n=3, 原子間: n=1の波動関数を構築する。C 原子: n=3, 原子間: n=1 +分極C 原子: n=2, 原子間: n=1 +分極+イオン

炭化水素系有機化合物 C 原子: n=3, 原子間: n=1 +分極+イオン
では、C-H 結合と C-C
Exact

Mn

Energy (a.u.)

426
2750
2902
1503
3790

-38.407 37-38.447 37-38.457 22-38.476 09-38.475 42-38.479ΔE(kcal/mol)
44.9
19.8
13.7
1.83
2.24

H-squareerror
4.71  10-1
2.00  10-1
2.01  10-1
4.06  10-1
1.80  10-1

結合がすべてであり、これらの効率的な記述を探求した。表 2 に、FATM-iExg 法による CH 分子の計算結果を示す。炭素原子の波動関数、原子間の関数、分極関数、イオン配置をテストした。
低オーダーの原子間の関数のみでは結合領域の記述が不足しているが、分極関数と特にイオン配置が効率的な記述に大きく寄与した。比較的少ない次元でも、絶対エネルギーで 2 kcal/mol 程度の誤差の解が得られた。当日は、C-C 結合も含め、詳細を発表する。
4. iExg 法による計算の加速
FATM-iExg 法は、2010 年 6 月にアイディアが提唱され、2013 年 4 月に東工大 TSUBAME グランドチャレンジ課題にて、高い並列パワーを駆使してプログラム開発と C2, N2, ベンゼン(C6H6) などのテスト計算を遂行した。表 3 に C2 分子の全電子の置換演算と iExg 法での置換演算を比較した(スピンから来る置換演算も考慮している)。前者は266!6!=33177600 個もの置換演算数を含み、我々が表 3. C2 分子の iExg 法による置換演算
交換する 交換の演 iExg 法の置換 全置換数に対電子数 算子の数
演算の数
する割合 (%)
01
384
0.0012
1
49
15744
0.047
2
353
1140480.34
3897
2891520.87
41201
3874561.17
51249
4028161.21
61250
4032001.22
全電子の置換演算

33177600

100.0

提唱したスレーター行列式の性質を利用する反対称化法(N3-5 の計算オーダー) [3] を利用しなけ
れば現実的な計算は難しい。一方、iExg 法では、隣り合う原子間で 3 電子程度の交換までで十分であり、その場合は電子交換の組み合わせが 897 通り、必要な置換演算は全体の僅か 0.87%で済む。ベンゼン、C60 など、分子がさらに拡大するほど飛躍的な計算の加速が見込まれる。
LSE 法の高い並列化効率はほぼ自明であるが[4]、TSUBAME でのベンゼンのテスト計算では超パラレルの結果が得られ、高い並列化性能が確認できた一方、単体でのチューニングの必要性が示唆された。TSUBAME での N2 分子とベンゼン分子の計算時間から計算オーダーを推察すると
O(N1.4)が得られ、大きな系で Order-N に近づくことが期待できる。現在、プログラムの改良と効率的な計算法の開発を進めている。当日はより大きな分子系の計算結果を発表する。
References: [1] H. Nakatsuji, J. Chem. Phys. 113, 2949 (2000). H. Nakatsuji, Phys. Rev. Lett. 93, 030403(2004). H. Nakatsuji, Phys. Rev. A 72, 062110 (2005). H. Nakatsuji, H. Nakashima, Y. Kurokawa, and A.
Ishikawa, Phys. Rev. Lett. 99, 240402 (2007). [2] H. Nakatsuji, Acc. Chem. Res. 45, 1480 (2012). 中辻 博,“最先端研究:「シュレーディンガー精度のシミュレーションを目指して」”, シミュレーション, 32,
39 (2013). [3] H. Nakashima and H. Nakatsuji, J. Chem. Phys. 139, 044112 (2013). [4] H. Nakashima, A.
Ishikawa, Y. Kurokawa, and H. Nakatsuji, SC12 conference in Salt Lake City, USA (2012).

4E08
Accurate and highly efficient approximations
approximations to
the twotwo-electron spinspin-orbit coupling terms(Institute for Molecular Science1)

Jakub Chalupský1, Takeshi Yanai1

Spin-orbit coupling (SOC), as an interaction of electronic spin with the magnetic
field caused by its movement relative to other charged particles of the system, isbelieved to be the most important spin-dependent relativistic effect. It gives rise,or significantly contributes, to many interesting physical and chemical
phenomena. Among those, we may mention for example intersystem crossing andspin-forbidden chemical reactions, phosphorescence, zero-field splitting, andelectron paramagnetic and nuclear magnetic resonance spectra. It is thus evident
that we may often need to include SOC in our quantum chemical calculations.
However, that is generally a non-trivial task, since the evaluation of two-electron
SOC integrals requires substantially, roughly by an order of magnitude, morecomputational effort than the evaluation of electron-electron repulsion integrals.
In this talk, we will present a new approximation to the two-electron SOC terms,
which we call flexible nuclear screening spin-orbit (FNSSO). This approximation
uses an effective one-electron SOC operator and accounts for the effect of twoelectron SOC by screening the nuclear charges involved in the corresponding oneelectron SOC integrals. Although this idea is not new (it is used for example ineffective nuclear charge and screened nuclear spin-orbit approximations), ourapproach takes full advantage of highly flexible scheme of the screening, wherebasically any pair of atomic orbitals has its own screening parameter. The
screening parameters come mainly from parameterization based on ab initio
atomic SOC calculations. Additionally, we use fairly simple physical model for
screening of multi-center interactions. As we will show on several practicalexamples (see Table 1 for brief summary of presented results), high flexibility ofthe screening used in FNSSO allows us to achieve mostly spectroscopic accuracy(error within 1 cm-1) for SOC matrix elements in molecules composed of lightatoms, and results with deviation of typically few wavenumbers for heavy-atom
containing systems.
In addition to its high accuracy and general applicability, FNSSO approach ishighly efficient, since only one-electron SOC integrals have to be evaluatedexplicitly. Moreover, it is also easy to implement, because only fairly simple
modifications of existing codes for one-electron SOC calculations are needed.

4E09
構造ゆらぎと化学反応に関する理論的研究	
 (阪大院・基礎工)	
  ○ 重田育照

Theoretical studies on chemical reactions with structural fluctuation(Graduate School of Engineering Science, Osaka University) ○Yasuteru Shigeta	  小さな孤立分子からナノマテリアルへと対象の広がりを見せる現代の分子科学は、より詳
細に分子の動的な側面を明らかにすることも一つの大きな研究の流れとなってきている。つまり、従来の静的な分子構造・物性・反応の独立した概念から脱却し、それら３つを複合した時空間の分子機能制御へと問題が複雑化している。その複雑な系の動的特性から基本原理を抽出する事が理論分子科学に与えられた大きな使命であり、その解析手段が量子化学的分
子動力学シミュレーションである。ナノスケールの領域での変化を記述する際には、小分子で重要であった、決定論的な「反応経路」が指標になるのではなくではなく、（非）平衡過程としての「動力学的効果」が重要な役割を果たすものと考え、これまで様々な研究を行ってきた。本講演は第５回分子科学会奨励賞記念講演として、これまでの研究を振り返りつつ、研究の将来像について述べる。以下では、代表的な幾つかの業績をまとめる。
（１）分子内包フラーレンの動的閉じ込め効果（核の温度揺らぎと電子物性変化）	  近年の合成技術の進展によって、籠状分子やメソポーラスやNaドープゼオライト等、ナノ空間を利用する様々な系が生み出されている。このようナノ空間に閉じ込められた分子には内部空間に於ける運動の自由度が存在する。その特徴は、多くの準安定構造とそれらをつなぐ反応経路が無数に存在する点である。我々は分子とその複合体の示す動的ゆらぎを計算科
学によって解明し、ダイナミクス-機能相関という新しい概念を構築したいと考えた。
	  そのような動的な閉じ込め効果が重要な系として、フラーレン内に水素分子とBe原子が内包されている系の理論計算を行った。特に、C60内でのBe原子の動的な振る舞いを調べるため、有効電荷の時間変化を半経験的QM/MM分子動力学法により解析し、温度と有効電荷や内部運動の関係を明らかにした。内包水素分子数がn=1（図1a）の時には電荷揺らぎの温度変化も小
さく、水素分子数がn=2（図1b）の場合には全ての温度領域において電荷揺らぎが非常に大きく、さらに温度変化によって、その分
布のピーク位置が変化することを明らかにした。また、その原因は形状に特徴的な運動と密接に関係すると言う、閉じ込め系において普遍的な現象を見いだした[1, 3]。	
  	
  	
  	
 
	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
 

図１：C60内でのBeのPESと有効電荷の動的相関	 (a) n=1, (b)n=2

（２）キュミュラント動力学理論の構築（核の量子揺らぎ・熱揺らぎの基礎理論）
	  分子はゼロ温度においても静止しておらず、零点振動をしている。これは最も基本的な量
子力学的効果である。しかし、動的な現象を解析するために現在広く行われている分子動力学計算では、核の運動は古典力学を用いて記述しているため、核の量子効果については範疇外であり、IRやラマン等の振動分光で捉えられる分子内振動すら近似的にしか取り扱う事は出来ない。そこで、我々は準量子キュミュラント動力学法を提唱した。この理論は、座標・運動量に加え、新たに量子揺らぎ（もしくは熱ゆらぎ）の変数を取り入れた、古典力学の自然な拡張になっており、今までの概念も包含している。例えば静的な理論として、量子論的な拡張変数を導入した	  ポテンシャル超曲面解析による量子トンネル経路の探索や、基準振動解析の確立、また動的な理論として、核酸塩基対の水素結合における軽水素と重水素置換体の動力学的安定性（量子動的同位体効果）を解析する方法論を確立した[2, 3]。最近は、古典統計力学的アンサンブルを評価する方法として、本手法を拡張している。
（ ３ ）自 由 エ ネ ル ギ ー 解 析 に よ る 複 雑 系 化 学 反 応 解 析（ 反 応 に お け る 揺 ら ぎ の 重 要 性 ）	  タンパク質や固体表面・界面での複雑で多段階の化学反応、ソフトマターや自己組織化複合体などの構造が柔軟な系における化学反応は、実験的にも非常に注目を集めており、反応過程の詳細において分子描像が求められている 。しかし、反応に寄与する自由度が多くなると、化学反応においても揺らぎの効果が無視出来なくなり、従来行われている静的なポテン
シャルエネルギー曲面解析を超えた動力学的な取り扱いが必要となる。我々は、メタダイナ
ミクス法に基づく自由エネルギー解析を様々な系に対して行い、複雑な化学反応の詳細な解析を行った[4,5]。例として、図2右側の2,4 Pentanediol-Tetheredは分子内ケテン-オレフィン[2+2]環化付加反応を起こし、99%の割合で図左側に示した一方の異性体が生成する事が実験によって報告されている。通常のポテンシャルエネルギー解析では、２つの異性体のエネルギー差は殆どなく、99%の選択性の理由は不明であった。メタダイナミクス法に基づく自由エネルギー面解析により、柔軟な構造を持つオレフィン部分の動力学的分布（エントロピーの効果）が、片方の選択的生成に寄与する事を明らかにし、反応における揺らぎの重要性を指摘した[5]。

	 図2 自由エネルギー解析と反応座標における相空間分布[参考文献] [1] Y. Shigeta, K. Takatsuka, J. Chem. Phys., 123 (Comm.), 131101 (2005). [2] Y. Shigeta,
H. Miyachi, K. Hirao, J. Chem. Phys., 125, 244102 (2006). [3] Y. Shigeta, Bull. Chem. Soc. Jpn. 82(Award Accounts), 1323 (2009). [4] K. Kamiya, Y. Shigeta, Biochim. Biophys. Acta., 1807, 1328-1335(2011). [5] K. Kamiya, T. Matsui, T. Sugimura, Y. Shigeta, J. Phys. Chem. A, 116, 1168 (2012).