1P-101

ORMAS-MCSCF に基づく相対論的多参照摂動法の開発(九大院理) 鈴木聡, 渡邉祥弘, 中野晴之

Multireference Perturbation Theory with Four-Component ORMAS-CI/MCSCF Reference Function(Kyushu univ.) Satoshi Suzuki, Yoshihiro Watanabe, Haruyuki Nakano【序】重原子を含む分子を記述するのに適した方法として、Dirac-Fock-Roothaan 方程式を出発とした多状態の多参照摂動法である相対論的 GMCQDPT[1-3]がある。これまでの研究ではいくつかの原子の励起スペクトルの計算を行い、全角運動量量子数の違いによる副準位のエネルギー差までよい精度で計算できることを確かめてきた。スピン多重度による制限がかかる非相対論の場合と異なり、相対論的な MCSCF ではすべての Slater Determinant を用いるため計算負荷は大きくなりがちである。また、励起状態を扱うために多くの状態を同時に扱えばその分計算負荷は増える。
そのため、大きな分子はもちろん原子の場合ですら、計算負荷が大きすぎ計算が実行できない場合があった。例えば、Gd 原子の基底状態(9Do2)を 4f,5d,6s からなる 26spinor を活性空間に用いた
CASSCF ではその次元は 500 万を超える。また、遷移金属の場合には(n+1)d 軌道まで活性空間に入れなければならない場合があり、この場合も CASSCF では取り扱えない。そこで、コンパクト
でありながら重要度の高い Determinant をとりこぼしのないように含まれる CI/MCSCF 波動関数を参照波動関数とすることで計算負荷を低減し、より大規模な系に対しても相対論的 GMCQDPT が適用できるようにすることを考える。
本研究では、Occupation restricted multiple active space (ORMAS) CI/MCSCF 波動関数[4]を相対論的波動関数に拡張し、参照波動関数に採用する。ORMAS-CI/MCSCF は活性空間を部分空間に分割し、各部分空間の占有電子数を指定することで重要な電子配置を効率よく生成することができる方法である。相対論的 GMCQDPT のプログラムに ORMAS-CI/MCSCF を実装し、ランタニド原子や遷移金属原子の励起状態を計算した。この予稿では Gd と Cu の例を示す。
【理論】GMCQDPT では、以下であらわされる有効ハミルトニアンを対角化することでエネルギ
ーを得る。

H

 E

(0 2)eff H I I H
 H I I H  
1


(0)
(0)
2
E  EI
E(0)  EI (0) I GCS I GCS    

MCSCF



ここで、相対論的ハミルトニアンの場合には、運動エネルギーが正の解のみを電子相関の計算に使用する、no-virtual-pair 近似を用いる。参照波動関数は任意の MCSCF あるいは CI 波動関数を用いることができる。
非相対論的な ORMAS では活性空間を任意の数の部分空間に分割し、それぞれの部分空間の占有電子数の最小値、最大値を指定する。この占有電子数に関する条件と、スピン多重度に由来する
α,β それぞれの電子数の条件を満たす Slater Determinant を CI/MCSCF の基底に用いる。これをスピンに関する制限を外すことで相対論的な波動関数に拡張する。スピノルで構成される活性空間
を部分空間に分割し、それぞれの占有電子数を指定する。そして、この条件を満たす Slater
Determinant 全てを CI/MCSCF の基底に用いる。
【 計 算 方 法 】 相 対 論 的 ハ ミ ル ト ニ ア ン と し て Dirac-Coulomb ハ ミ ル ト ニ ア ン を 用 い 、Dirac-Fock-Roothaan SCF は DIRAC11[5] により行った。
Gd 原子の基底状態 9Do2 である。9Doj(j=3,4,5,6)への励起エネルギーを状態平均 MCSCF および

GMCQDPT により求める。これらの電子状態は(4f)7(5d)1(6s)2 に対応する。そこで、Active Space を
f 軌道、d 軌道、s 軌道からなる三つの subspace に分割し、それぞれに丁度 7,5,2 個ずつ電子が入る
ような全ての Determinant を用いて計算を行うことにする。このようにすると、MCSCF の次元は
CASSCF を用いた場合の 1/100 未満の 34320 次元になり、容易に計算可能になる。
基底関数は Dyallの基底を(8s7p5d4f)に短縮して用いた。内殻の 36 電子は電子相関の計算には用いず、5s 軌道以降を相関させた。
Cu の基底状態は(3d)10(4s)1 に対応する 2S1/2 と、(3d)9(4s)2 に対応する 2Dj(j=5/2,3/2)を状態平均 MCSCF
および GMCQDPT により求めた。3d,4s,4p,4d からなる 28spinor を活性空間とし、3d 軌道に 7～9電子、4s 軌道に 0～2 電子、4p,4d からなる subspace に 0～2 電子が入るような全ての Determinantを用いた。
これは 3d,4s からなる 12spinor を活性空間に用いた CASSCF からの 2 電子励起を含む、
MRCISD のような波動関数である。比較のために 12spinor を活性空間に用いた CASSCF による計算も行った。基底関数は古賀らの基底関数[6]を短縮し、分極関数として f 関数を追加し、(6s5p3d2f)とした。この計算ではすべての電子を相関させた。
【計算結果】Gd の計算結果を表 1 に示す。CASSCF に較べてずっとコンパクトであるにも関わらず、9Dj 状態の記述に必要な最低限の Determinant は適切に含まれているため、実験値とよく一致している。より高い励起状態は 11Foj は 6s→5d 励起を含んでいるので、対応する Determinant を含めることで取り扱える。一方、表 2 に示したように Cu については 12spinor を活性空間とした
CASSCF では 2D5/2 が最低エネルギーとなり定性的に間違った結果を与える。摂動により振る舞いは改善するが実験値を再現しない。一方、ORMAS を用いて MRCISD 様の参照波動関数を用いた
場合、MCSCF の段階で 2S1/2 が正しく基底状態となっている。また摂動により実験値を再現する。
このように、重要な電子配置を的確に選択することによって、コンパクトな波動関数でも精度よく電子状態を記述できる。
MCSCF

GMCQDPTo

D2

0.00000.00000.00009Do30.02370.03010.02672

S1/2

2.87430.00000.00000.00000.0000

9

Expl.

Term

CAS(11e.,12s.)
CASSCF MCQDPTORMAS3(11e.,28s.)
MCSCF GMCQDPT

Expl.

9

o

D4

0.06090.06490.06612

D5/2

0.00001.91220.22551.3294

1.38899oD50.11940.15470.12392

0.27852.17520.50371.61481.64229Do60.21400.29870.2131

D3/2

表 1. Gd の 9Dj 状態

表 2. Cu の d→s 励起【参考文献】[1]Miyajima, M.; Watanabe, Y.; Nakano, H.; J. Chem. Phys. 2006, 124, 044101.
[2]Ebisuzaki,R.; Watanabe,Y.; Nakano,H.; Chem. Phys. Lett.2007,442, 164.
[3] Ebisuzaki,R.; Watanabe,Y.; Kawashima,Y.; Nakano, H.; J. Chem. Theory Comput. 2011, 7, 998.
[4]Ivanic,J; J. Chem. Phys. 2003, 119, 9364.
[5]DIRAC, Release DIRAC11 (2011) Bast R., Jensen H.J. Aa., Saue T., Visscher L.
[6] Koga, T; Tatewaki, H; Matsuoka, O;. J.Chem.Phys. 2001, 115, 3561/

1P102有機 EL 錯体における燐光過程の理論的考察（阪府大院理 1, RIMED2, JNC Co.3）小関史朗 1,2, 鎌田尚也 1, 鍵田侑希 1,
麻田俊雄 1,2, 松下武司 2,3Theoretical Investigations on OLED Phosphorescent Complexes（OPU1, RIMED2, JNC Co.3）Shiro Koseki1,2, Na-oya Kamata1, Yuhki Kagita1,
Toshio Asada1,2, Takeshi Matsusita2,3【序】
液晶ディスプレイに代わる次世代ディスプレイとして注目を集めてきた有機 EL ディスプレイは，価格，寿命および効率の問題により期待されたほど普及していない．
今後の普及のためには，その技術開発には何らかの大きな変革が必要とされている．
本研究シリーズでは，燐光を発する白金およびイリジウム錯体に着目し，燐光波長への配位子の効果および配位子への置換基導入の効果を理論的に研究し，変革をもたらす新規発光材料の提案を試みる．
【計算方法】
基底状態および最低三重項状態の幾何学的構造を密度汎関数法（B3LYP/SBKJC+p）により最適化し，適切で実行可能な活性空間を用いた MCSCF 法により分子軌道を求める．さらに，動的相関を考慮した Second-order configuration interaction 法により波動関数を構築し，それらを用いてスピン軌道相互作用行列を構築・対角化し，電子状態（Spin-mixed (SM) states）の波動関数を求めた．なお，ECP 基底関数を用いたため，スピン軌道相互作用の見積もりには Breit-Pauli Hamiltonian の一電子近似を用いた．これらの電子状態間の遷移モーメントを計算することにより，燐光スペクトルの波長と強度を予測した．
【結果と考察】
N
N
●白金錯体
Pt
親分子として Pt(thpy)2 を用い，配位子の置換，配位子への置SS換基の導入に伴う燐光スペクトルの波長および強度の変化を見
cis-Pt(thpy)2積もった．特に，塗布型プロセスによるデバイス作製を目的とした溶解度の高い Thompson 型錯体(C^N)Pt(O^O)に着目し，より鮮明な赤色の燐光を発する錯体の設計を目指す．本研究では，Pt(thpy)2 
thpyPt(acac)  thpyPt(bpp)  thpyPt(bdmp)のシリーズについて理論計算を行い，下記の項目が燐光波長に与える影響について解析した：
１．acac 配位子の効果，
２．acac 配位子への MeO 基導入の効果，
３．acac 配位子への MeO 基導入位置，

４．thpy 配位子へのベンゼン環導入の効果．
５．ppy 配位子の効果．
これらの結果は論文において報告する準備を進めている．詳細は当日報告する．
●イリジウム錯体
白金錯体と同様に，親分子として Ir(ppy)3 を用い，配位子のN置換，配位子への置換基の導入に伴う燐光スペクトルの波長おIr
よび強度の変化を見積もった．塗布型プロセスによるデバイス
NN
作製を念頭におき，O^O 配位子導入の影響を考察した．イリジ
ウム錯体では，特により鮮明な青色燐光を発する錯体の提案・合成が期待されており，多々研究が進められている．現在まで
fac-Ir(ppy)3に最も有用であるとされているのは，フッ化物を導入した
Ir(dfppy)2(pic)，いわゆる FIrpic であり，これ以上に有用な錯体の提案が強く望まれている．本研究では，Ir(ppy)3  Ir(ppy)2(dfppy)  Ir(ppy)(dfppy)2  Ir(dfppy)3 Ir(dfppy)2(pic)  Ir(dfppy)2(acac)のシリーズについて理論計算を行い，下記の項目が燐
光波長に与える影響について解析した：
１．ppy 配位子へのフッ素基導入の効果，
２．ppy 配位子へのフッ素基導入位置の効果，
３．pic 配位子導入の効果，
４．acac 配位子導入の効果．
なお，Ir(ppy)2(dfppy)や Ir(dfppy)2(pic)のような Ir(C^N)2(C’^N’)あるいは Ir(C^N)2(N^O)
では４つの構造異性体（HNT, HCT, HC および HT），Ir(dfppy)2(acac)のようなIr(C^N)2(O^O)では３つの構造異性体（HNT, HCT および HC）が存在する．光学異性体も存在するが，本研究で用いる分子軌道計算では相違が得られないので考察から除いた．我々の計算によれば，どの構造異性体を用いるかによって燐光波長が大きく変化する．それゆえ，これらの構造異性体の異性化反応が燐光波長に与える可能性についてもモデル錯体を用いて考察したので報告する．これらの結果をもとに適切な青色
燐光材料分子の提案を行っていく．これまでの結果に関して，現在，論文発表の準備中である．
【参考文献】
分子科学討論会要旨（1P117(2010), 3P111(2010), 2P111(2011), 4P109(2011)）．
理論化学討論会要旨（2P24(2012)）．
T. Matsushita, T. Asada, S. Koseki, J. Phys. Chem. A, 2006, 110, 13295.
T. Matsushita, T. Asada, S. Koseki, J. Phys. Chem. C, 2007, 111, 6897.

1P-103
時間依存密度汎関数法に基づく sum-over-states 法による
p-ニトロアニリンの電子遷移特性と第一超分極率についての理論的研究（阪大院基礎工）村田裕介、岸本真悟、岸亮平、中野雅由

Theoretical study on the electronic transition properties andthe first hyperpolarizability of p-nitroaniline by the sum-over-states method
based on the time-dependent density functional theory(Graduate School of Engineering Science, Osaka Univ.) Yusuke Murata, Shingo Kishimoto,
Ryohei Kishi, Masayoshi Nakano【序】時間依存摂動論に基づく sum-over-states (SOS)法では、分子系の静的第一超分極率の z-軸成分βzzzは以下のように表される。

µ z (µ z ! µ z " )µ z
! zzz = 3## 1n nm 11 nm m1
En1Em1n"1 m"1(1)zz
ここで En1 , µ nm
, !µ nz = µ nn
! µ11z はそれぞれ状態 n の励起エネルギー、状態 n, m 間の遷移モーメント、

状態 n と基底状態(1)の双極子モーメント差を表す。SOS 法はこれら電子状態間の遷移特性（電子遷移特性）と非線形光学応答とを直接関係づけることができるため、二次非線形光学物質の構造-特性相関について有益な情報を与える。実際、有限場(FF)法や応答理論などによるβの計算結果の解析において、式(1)で主な寄与を与える 2 状態のみを考慮する近似（2 状態近似）などが広く適用されている。一方、SOS 法の結果は、一般に励起状態計算による電子遷移特性の記述の正確さや考慮する状態数 N に依存するため、使用に際してはこれらの依存性を検討する必要がある。近年、励起状態計算には、低い計算コストで比較的高精度な結果が得られる時間依存密度汎関数 (TDDFT) 法が広く用いられている。
しかしながら、これら電子遷移特性の記述のバランスが顕著に結果に反映する非線形光学応答量のSOS 計算への適用性や、その汎関数依存性についての検討は殆どなされていない。そこで本研究では、
対象分子として p-ニトロアニリンを取り上げ、高精度 ab initio MO 励起状態計算法と TDDFT 法の電子遷移特性と静的β の長軸方向成分(βzzz)の計算結果を比較し、これらの物性量に関する SOS 法の適用性を議論する。
【計算】p-ニトロアニリンの構造最適化は B3LYP 法により行った。本研究では、Tamm-Dancoff 近似
における TDDFT (TDA)法を用い、汎関数依存性を検討する。汎関数は、B3LYP 法の他に、長距離補正(LC)-BLYP 法を適用する。LC 法の領域分割パラメータは、DFT 法のエネルギー線形条件と Janak の定理に基づく IP-tuning 法[1]により最適化した 0.297 bohr-1 を用いた。次に、双極子モーメント差については以下の二種類の定義を検討する。
軌道緩和なし（unrelaxed）の双極子モーメント差：!µ nUD = ! µ! n (r)dr " ! µ!1 (r)dr(2)

軌道緩和あり(relaxed)の双極子モーメント差：!µ nRD = !

"En1 ( F )"F F=0(3)無電場下での状態 n の期待値として定義する unrelaxed 双極子モーメント(UD)は、分極の無摂動系の固有関数による展開に基づく SOS 法においてよく用いられる。一方、relaxed 双極子モーメント(RD)は、
印加電場による軌道緩和の効果を取り込むため、特定の励起状態の双極子モーメントの記述に優れると期待されるが、様々な汎関数を用いた場合の DFT 計算における信頼性は未解明である。参照として、
ab initio MO 励起状態計算法である CIS 法、EOM-CCSD 法の結果と比較した。高精度電子相関を取り込める EOM-CCSD 法では、RD の場合のみを検討した。また CCSD 法に基づく有限場(FF)法の結果とも比較した。全ての計算には、6-31G**基底関数系を用い、GAUSSIAN09、GAMESS および NWChem 6.0
プログラムパッケージにより行った。
【結果と考察】まず、長軸方向に基底状態から零でない遷移モーメントを持つ最初の励起状態 k についての 2 状態近似の結果を表 1 に示す。B3LYP で過小評価していた励起エネルギーは長距離のHartree-Fock 交換項を考慮することにより改善する結果を得た。遷移モーメントについては汎関数依存性が小さく、B3LYP 法と LC-BLYP 法はほぼ同じ値を与える。一方、双極子モーメント差（UD）は、
大きな汎関数依存性を示した。TDA 法の UD は EOM-CCSD の結果(RD)を著しく過大評価し、大きさは CIS (UD) [0.73] < EOM-CCSD (RD) [1.0] < LC-BLYP (UD) [1.45] < B3LYP (UD) [2.02]の順となった（[ ]内は比率を表す）。これよりβ値も同様に過大評価されることがわかった。これに対して RD の結果は逆に過小評価し、大きさは CIS(RD) [0.56] < B3LYP (RD) [0.63] < LC-BLYP (RD) [0.75] < EOM-CCSD(RD) [1.0] の順となり、UD に比べて汎関数依存性が小さく、β値は EOM-CCSD の値に近くなった。
	  次に SOS 法で考慮する励起状態数に対するβの依存性について、式(1)において n = m の対角項のみ考慮する近似を適用した結果を図 1 に示す。どの手法においても UD を用いた結果に比べ RD を用いた結果は、考慮する状態数や汎関数に対する依存性が小さくなることが明らかになった。RD を用い、20 状態を考慮した場合のβの結果は、CIS (592 a.u.) < EOM-CCSD (852 a.u.) < LC-BLYP (910 a.u.) <B3LYP (956 a.u.) であり、LC-BLYP 値が EOM-CCSD 値に最も近くなった。以上の結果より、p-ニトロ
アニリンにおいて SOS 法の対角近似を用いる場合、RD の定義を用いた LC-BLYP 汎関数の結果が電子遷移特性とβ値をともに最も良く再現することが明らかになった。発表では、基底状態と励起状態の差密度および遷移密度を解析し、UD における汎関数依存性の原因についても検討する予定である。
表 1. 電子遷移特性および 2 状態近似での第一超分極率β| µ1k | [D]

EOM-CCSD

4.889

4.732

CIS

5.454

4.900

4.757

5.261

4.280

5.266

–

–

TDA/
LC-BLYP
TDA/
B3LYP
CCSD FF 法aΔµk [D]β [a.u.]

–
6.865
5.000
3.855
9.970
5.137
13.90
4.319

–
870
546
421
1650
885
2847
839

–

817

EOM-CCSD (RD)LC-BLYP (RD)LC-BLYP (UD)B3LYP (RD)

a

3000

! [a.u]Ek1 [eV]a

4000

上段: UD,下段: RD【参考文献】 [1] T. Stein et al., J. Am. Chem. Soc., 131, 2818 (2009).

B3LYP (UD)CIS (RD)CIS (UD)

2000
1000
0

図 1.

0

5

10
N

15

SOS 法によるβの状態数(N)依存性20

1P104量子波束法を取り込んだ非断熱遷移を扱う第一原理分子動力学法の開発(北大院・総合化学 1, 北大院・理 2) 小林 孝生 1, 野呂 武司 2, 武次 徹也 2

AIMD with nonadiabatic transition estimated by the wavepacket approach(Hokkaido Univ.) Takao Kobayashi, Takeshi Noro, Tetsuya Taketsugu【研究背景】
近年、光反応のメカニズムを理解するために励起状態ダイナミクス計算が行われるようになってきたが、複数の電子状態が近接し断熱近似が成り立たない領域では、断熱状態間の非断熱遷移を考慮しなくてはならない。化学的に重要な過程には非断熱遷移を伴うものが多く、第一原理分子動力学法(ab initio molecular dynamics(MD))の計算において非断熱遷移を取り入れる効率的かつ実用的な方法が求められている。
ab initio MD 計算に非断熱遷移を取り入れる方
法としてもっともよく用いられているのは Tullyの状態遷移アルゴリズム[1]である。この方法の問
題点の一つは、MD の各ステップで状態間の遷移を見積もるためにコストの大きな「非断熱結合」の計算を要する点である。そこで本研究では、
非断熱結合の計算回数を減らすことで計算コストの軽減を図ることを目的として、非断熱遷移 図 1 Zhu-Nakamura 理論と一次元の切り出しに関する Zhu-Nakamura 理論と量子波束法を組み合わせた方法を考える。Zhu-Nakamura
理論では状態間の非断熱遷移はポテンシャ
ル曲面が最も近接した点で起こると考え、その近傍のポテンシャルの形状によって求められる遷移確率により非断熱遷移の有無を判定する。この際、非断熱遷移の方向を決定

図 2 Zhu-Nakamura 理論で判定可能な例

するために最近接点において非断熱結合の計算が必要となる。つまり、計算コストが高い非断熱結合はポテンシャル曲面の最近接点においてのみ必要となり、計算時間の大幅な軽減が期待される。Zhu-Nakamura 理論は 1次元のポテンシャル曲線に対して定義された理論であるので、多次元である ab initio
MD にそのまま組み込むことはできない。そこ

図 3 Zhu-Nakamura 理論で判定不可能な例

で、MD でポテンシャルの最近接点に到達したら非断熱結合ベクトルを計算し、トラジェク
トリとは別にそのベクトルの方向に 1 次元のポテンシャル曲線を計算し、それに対して
Zhu-Nakamura 理論を適用する[2]。

しかし、Zhu-Nakamura 理
論の公式は Landau-Zener 型と Nonadiabatic Tunneling 型の 2 種類のポテンシャル曲
線にしか適用できない。そこ
で、Zhu-Nakamura 公式が適用できない場合に対して、補完的に量子波束法 [3] を適用することにした。この量子波

図 4 量子波束法の概要

束法は、前述の 1 次元のポテンシャル曲線上でガウス型の波束を走らせ、各状態への波束の分岐比から非断熱遷移の確率を求める方法である。ここで実際に波束が走るのは電子状態計算で求めた断熱ポテンシャルではなく、状態同士が交差する透熱ポテンシャルである。あらかじめ 2 状態の透熱ポテンシャル(Vii、Vjj)をパラメータが含まれる関数で表わし、Vii を対角成分に、Vij を非対角成分に持つ行列を対角化した時に、元の断熱ポテンシャルになるようにパラメータを決めたものである。この透熱ポテンシャル Vii をもとに、ハミルトニアンによる時間発展演算子を波束に対して演算する。次に波束の状態遷移についての演算として、カップ
リング Vij に関する演算子を演算する。これを繰り返して波束を時間発展させていく。

【結果と考察】この量子波束法を ab initio MD に組み合わせて、シス体-トランス体間の光異性化反応を起こすアゾベンゼンをテスト分子として用いた。第一励起状態(S1)を初期状態とし、基底状態(S0)への非断熱遷移を評価する。いくつかの異なる初期構造の
うち Zhu-Nakamura 理論が適用可能ない
くつかの例について、遷移確率の結果をまとめたものが図 5 である。横軸に
Zhu-Nakamura 理論、縦軸に量子波束法に 図 5 Zhu-Nakamura と量子波束法の遷移確率の比較
よる遷移確率をプロットした。遷移確率が大きい領域では両者に良い一致が見られるが、遷移確率の小さい領域では多少の差が見られる。より詳細な議論は当日報告する。

【参考文献】[1] J. C. Tully, J. Chem. Phys. 93, 1061 (1990).
[2] S. Nanbu, T. Ishida, and H. Nakamura, Chemical Science 1, 663 (2010).
[3] X. Chen and V. S. Batista, J. Chem. Phys. 125, 124313 (2006).

1P105F-CN Bond Cleavage by Silyl Iron Complex. Electronic Effect of the Group or Atom
Attached to Cyano Group. A DFT Theoretical Study(Grad. Sch. Inf. Sci., Nagoya Univ.1, Fac. Sci., Assiut Univ.2)
AbdelRahman A. Dahy1,2, Nobuaki Koga1
1.Introduction:
Cleavage of strong bonds in some ligands will give new species which are important in organic
syntheses. There are a number of experimental and theoretical works focusing on cleavage of these
unreactive bonds, especially E−CN (E = C, N, O) bonds in the presence of iron or molybdenum silyl
complexes in cases of E = N and O by Nakazawa et al.1-4 The idea behind using transition metals is that themetal attaching to the ligand through its special part will affect the electronic properties of the strong bondand that the silyl group will facilitate this cleavage. As a matter of fact, in the case of CH3−CN, a new
species which contains CN and CH3 ligands are formed1 while in the cases of cyanamides and cyanates,
Et3SiCN is formed, to demonstrate the E–CN bond cleavage.2,3 The N−CN bond is difficult to be broken;
because, for example, in the case of H2NCN, the lone pair of electrons on amino N atom strengthens thisbond (because of its tautomerization to HN═C═NH), and it is even stronger than CH3−CN bond.2,4 Nsilylated η2-amidino iron complex was isolated and identified by x-ray crystallography in the case of
cyanamides and its heating lead also to the formation of Et3SiCN.2 Accordingly, Nakazawa el al. suggested
the mechanism shown in Scheme in the cases of cyanamides as well as cyanates. It was found that the low
activation energy for the silyl group migration from Fe to NCN is attributed to the hypervalent character of
Si atom and the charges on the Si and N atoms.
In this theoretical study, continuing to our study for the cleavage of E−CN bond in the presence ofsilyl iron complex; (i) the cleavage of F−CN will be studied according to Scheme to show if they are
similar to those in the cases of E = C, N, and O; and (ii) the effect of the group attached to CN on thereaction mechanism will be studied in detail.
2. Computational methods:
Density functional theory (DFT) calculations were performed to verify the feasibility of the reaction
mechanism suggested in Scheme. The optimized structures of intermediates and transition states and the
profiles for potential energy corrected with zero-point energy (ZPE) and the Gibbs free energy for thestudied mechanism are shown in Figure. Unless explicitly stated otherwise, the ZPE-corrected potential
energies are used in the discussion.
3. Results and Discussion:
As shown in Figure, searching for a minimum for the complex of FCN with complex 1, leads to theend-on (η1-) complex, 2d. 2d is more stable than 1 and FCN by 16.9 kcal/mol. This is because the strong
coordination of FCN to 1 via the lone pair of electrons on the NCN atom. The end-on complex, 2d,isomerizes to the side-on (η2-) complex, 3d, via the TS2d-3d. The activation energy for this step is 13.2
kcal/mol. The vibrational mode with an imaginary frequency at TS2d-3d is the Fe····CCN bond distancestretching. Accordingly, this bond distance shortens during the reaction from 3.054 Å in 2d to 2.745 Å in
TS2d-3d, and the Fe−CCN bond forms in 3d with a bond length of 1.882 Å. The next step is the migration of
the silyl group from Fe to NCN which transforms 3d to 4d via TS3d-4d. The activation energy of 0.6 kcal/mol
for this step is relatively small. The vibrational mode with an imaginary frequency at TS3d-4d is the N····Sibond distance stretching. Accordingly, during the reaction this bond distance shortens from 2.695 Å in 3d
to 2.329 Å in TS3d-4d and then the N−Si bond forms in 4d with a bond distance of 1.746 Å. This small
activation energy is because of a large exothermicity of −32.6 kcal/mol, the hypervalent character of Si
atom, and the strong electrostatic interaction between the positively charged Si and negatively charged N
atoms as stated previously. The Mülliken charges on the Si and N atoms in 3d are +0.923 and −0.296,
respectively. Therefore, the Si atom may be strongly attracted by the N atom to form a strong bond. In thenext step, it is expected that 5d having the Fe−C−F three-membered ring, similar to 5's in the cases of
H2NCN, Me2NCN and MeOCN, is formed. In fact, 5d does not form from 4d but instead 4d is directlytransformed to 6d via TS4d-6d. The activation energy for this step is 15.9 kcal/mol and it is the ratedetermining step. The vibrational mode with an imaginary frequency is the combination of cleavage ofFe−N bond and formation of Fe−F bond. Accordingly, the Fe−N bond is elongated from 2.123 Å in 4d to
2.890 Å in TS4d-6d and then to 3.006 Å in 6d. On the other hand, the Fe−F bond shortens from 3.082 Å in

4d to 2.773 Å in TS4d-6d and eventually to 1.854 Å in 6d. This reaction coordinate at TS4d-6d is similar to
those at the TSs for 4→5 with E=C, N, and O, to suggest that neither 5d nor TS5d-6d exists. This is ascribedto the large stability of 6d; it is 23.5 kcal/mol more stable than 4d in contrast to the reactions with E=C, N,
and O in which 6a and 6c are more stable than 4a and 4c by only 5.9 and 7.9 kcal/mol, respectively, and 6band 6b’ are less stable than 4b and 4b’ by 5.9 and 10.1 kcal/mol, respectively.
In conclusion, the reaction mechanism with E = F is similar to those in cases of E = C, N, and O. The
activation energies for the rate-determining steps are 28.9, 29.9, 28.0, 19.1, and 15.9 kcal/mol in the casesof ERn = CH3, H2N, Me2N, MeO, and F, respectively. Those for MeO and F are much lower than theothers; while in the rate-determining step with ERn = CH3, H2N, and Me2N the structures of ERn groups
must change drastically, the lone pair of electrons on the O and F atoms can easily interact with the Fe
atom.4 The large exothermicity in the case of ERn = F makes a contribution to the lower activation energy.
Scheme. Proposed Scheme for E-CN bond cleavage in the presence of iron silyl complex.

[Fe]
ERn

ERn

3x

2x

N
C

[Fe]

C

C
ERn

SiMe3

C

[Fe]

[Fe]

N C ERn

N

N

N
[Fe]

SiMe3

SiMe3

SiMe3

SiMe3

ERn

4x

6x

5x

[Fe] = CpFe(Co)ERn = CH3 (x=a), NH2 (x=b), NMe2 (x=b'), OMe (x=c), and F (x=d)

8

30

2
1.

2d

63

TS2d-3d(b)TS2d-3d
3.7
3.4]

1.179

2.
44
1.231
25

29

TS3d-4dTS3d-4d
24.6
15.4]3d
25.2
16.7]TS4d-6d
41.9
35.9]

4d

9
1.76

1.771

1 .9

2

2. 7

73
35

6

1.
34

3
1.

27

74

1.
22 246
1.4

8
1.

1.

90
2.8

23
2.1

1.8
43

0.0
[0.0] 2d
16.9
10.0]

1.260

8
1.85

3d

1+ FCN

54
1.8

1.219

2.
3

29
1.3

6

2.4
12

66

95

45
2.7
1.1

2. 6

1.894
1.273 1.159

2
1.
90

1.274

1.162
1.1

0

91
2.0

1 + FCN

1.
93

06
2 .1

1.9
35

22

7
2.16

1
2.

1.882

H

2.3
55

F

2.

O

34

C N

Fe Si

9

(a)

57.8
52.5]6d
81.3
76.2]4dTS4d-6d

6d

Figure. The (a) optimized structures and (b) the energy profile for the cleavage of F−CN bond in thepresence of iron silyl complex. All bonds are in Ǻ. All energies are ZPE-corrected and relative to
CpFe(CO)SiMe3), 1, and F−CN. Energies in parenthesis are free energies relative to 1 and F−CN.
References:(1) Nakazawa, H.; Kawasaki, T.; Miyoshi, K.; Suresh, C. H.; Koga, N. Organometallics, 2004, 23, 117.
(2) Fukumoto, K.; Oya, T.; Itazaki, M.; Nakazawa, H. J. Am. Chem. Soc., 2009, 131, 38.
(3) Fukumoto, K.; Dahy, A. A.; Oya, T.; Hayasaka, K.; Itazaki, M.; Koga, N.; Nakazawa, H.
Organometallics, 2012, 31, 787.
(4) Dahy, A. A.; Koga, N.; Nakazawa, H. Organometallics, 2012, 31, 3995.

1P106
水素ハイドレートにおけるケージ構造の占有率依存性（明治大院・理工）伊東 雄，深澤 倫子

Dependence of cage occupancy on cage structure of hydrogen hydrate（Meiji Univ.）Yu Ito，Tomoko Ikeda-Fukazawa【序論】
クラスレートハイドレート（以降、ハイドレートとよぶ）は、水分子によって構成されるかご状の格子（ケージ）中にゲスト分子を取り込むことで安定化する包接水和物結晶である。ハイド
レートはその包接特性により、安定かつ高密度に気体分子を貯蔵可能であるため、エネルギー資
源となる気体分子の輸送媒体や貯蔵媒体として期待されている。近年特に、水素をゲストとする
水素ハイドレートが、クリーンエネルギーである水素の貯蔵媒体として注目される。
水素ハイドレートのホスト格子は、II 型構造と呼ばれる立方晶系（空間群 Fd3m）をとる。この II 型構造の単位格子は、136 個の水分子から成り、8 個の 16 面体（以降、L ケージとよぶ）と16 個の 12 面体（以降、S ケージとよぶ）によって構成される。ハイドレートは非化学量論物質であるため、ケージ中のゲスト分子の含有率（ケージ占有率）は、温度・圧力等の熱力学条件に依存して変化する。水素ハイドレートの場合、温度 160 K、圧力 200 MPa の条件で、最大のケージ占有率（L ケージの占有率（L）= 4、S ケージの占有率（S）= 1）になると報告されている[1]。
水素ハイドレートの実用化を実現する上で、現存の水素貯蔵媒体の水素重量比を超える高いケー
ジ占有率を有するハイドレートの生成・保管法を確立することが重要となるが、最大ケージ占有率の達成条件や熱力学条件の変化に伴う占有率の変化の過程には不明な点が多く残されている。
さらに、包接するゲスト分子の運動状態や占有率の変化に伴ってホスト格子に歪みが発生することが分かっているが[2]、複占有に伴うケージ歪みおよび歪みに起因した構造安定性の変化については明らかではない。
本研究では、ケージ間の分子移動の誘因となるケージ複占有のメカニズムを解明することを目的として、水素ハイドレートの分子動力学計算を行った。この結果を基に、ケージ占有率の変化に伴うホスト格子の構造歪みを解析した。
【実験】
本研究では、水分子 1088 個から成る 2×2×2 の 8 ユニットセルを格子とする計算系を用いた。ケージ占有率については、S を 1 に固定し、L を 1-4 の範囲で 1 刻みに変化させた。計算条件はNTP アンサンブルとし、圧力は 220 MPa、温度は 10K とした。
分子動力学計算には、プログラム MXDORTO[3]を用いた。原子間相互作用の計算には、KAWAMURA ポテンシャルモデル[4]を使用した。また、計算ステップは 0.4 fs/step とした。

【結果と考察】
ケージを構成する水分子の密度分布の変化を解析した結果、L

1.4

(a)

θL=1
θL=2
θL=3
θL=4

1.2

ケージ、S ケージ共に、L ケージの占有率の増加に伴ってケージが膨張し、構造歪みが大きくなる
ことが明らかになった（Fig. 1）。
L ケージ中にゲストとして包接さ

1.0
0.8
0.6

れる水素分子は、複占有によりケージ中心から外側に外れて分布する。このことから、L ケージの

0.4
0.2

構造歪みについては、包接する水
素分子と水分子との反発相互作

0.0
4.0

用によって引き起こされたもの

4.2

4.4

4.6

4.8

5.0

5.2

Distance from L-cage center[Å]であると考えられる。一方、1 個の水素分子のみを包接する S ケージの膨張については、L ケージの膨張に伴う水分子の配置の変化に起因したものであると考えら

1.4
1.2

θL=1
θL=2
θL=3
θL=4

(b)

1.0

れる。
さらに、ケージを構成する水分子の熱振動の変化を解析した結

0.8
0.6

果、L ケージ、S ケージ共に、L
ケージの占有率の増加に伴い、熱振動の振幅が減少することが分

0.4
0.2

かった。高占有率の状態では格子の構造歪みが増大しているにも関わらず熱振動の振幅が減少することから、複占有のハイドレー
トに発生する格子歪みには構造安

0.0
3.2

3.4

3.6

3.8

4.0

4.2

4.4

Distance from S-cage center[Å]Fig.1 (a)L ケージ、(b)S ケージを構成する O 原子の密度分布。

定化の効果があると考えられる。
以上の結果を基に、占有率の増加に伴うケージの構造歪みのメカニズムについて議論する。

[1] K. A. Lokshin, Y. Zhao, D. He, W. L. Mao, H. K. Mao, R. J. Hemley, M. V. Lobanov, M. Greenblatt:, Phys. Rev. Lett., 93,(2004) 125503[2] T. Ikeda-Fukazawa, Y. Yamaguchi, K. Nagashima, and K. Kawamura, Journal of Chemical Physics 129 (2008), 224506.
[3] K. Kawamura, MXDORTO, Japan Chemistry Program Exchange, #029.
[4] N. Kumagai, K. Kawamura, and T. Yokokawa, Mol. Simul. 12 (1994) 177.

1P107
二重指数関数型公式による２原子分子の電子密度の数値積分（三重大院工）三谷 昌輝・吉岡 泰規
Numerical Integration of Diatomic Molecular Electron Density with
Double Exponential Formula(Mie Univ.) Masaki Mitani, Yasunori Yoshioka【序】密度汎関数計算では交換・相関エネルギーをエネルギー汎関数の数値積分で求めるため、計算結果は数値積分の精度に依存する。通常、多中心の分子積分は１中心の原子積分の和に分解され、原子積分は原子核を中心とする極座標を用いて動径変数 r と角度変数(, )に分離し、動径積分と角度積分を数値積分により近似計算する。これまでに、数値積分を行うための種々の積分グリッドが提案されているが、標準的な積分グリッドに対して幾つかの問題点も報告されている。
本研究では、動径関数の数値積分に対して二重指数関数型公式を適用し、その有効性を検証した。
【方法】これまでの研究では、動径関数の数値積分に対して、主に(1)式の第２種 Gauss-Chebyshev
公式や(2)式の Euler-Maclaurin 公式が用いられている。これらの求積公式を動径積分に適用する際に、マッピング変換 r = t(x)により動径変数 r を求積変数 x に変数変換するが、この変換には任意のパラメータが含まれる。マッピングパラメータを変化させることによりグリッド点の分布が変化するため、積分結果はグリッド点の数との値に依存する。

 1  x 
1

1

2 12

f  x  dx n

 i 

 sin  n  1  f  x 
n 1
2

(1)

i

i 1

n

1

 f  x  dx  n  1  f  x 
1

0

i 1

(2)

i

積分グリッドは、用いる求積公式とマッピング変換により定まる。代表的な動径グリッドとして、(1)式の求積公式と(3)式の変数変換を組合せた(4)式の TA グリッド[1]や(2)式の求積公式と(5)式の変数変換を組合せた(6)式の MK グリッド[2]が提案されている。
t  x  

1  x 

0.6

ln 2

1 x 
ln 

 2 

(3)

12
12
3

1  xi  3  1  xi 
1.8  1  xi     n2
2  1  xi 
F
r
r
d
r

1

x
ln

0.6
ln 
  i   1  x 1 2  2   F  ri 
12
0  
2
 ln 2  n  1 i 1


1

x





i
i




t  x    ln 1  x3 23
 F  r  r dr  3
0

(4)
(5)

1 n xi ln 1  xi  1  x3 F  ri n  1 i 1
i
2

2

3

(6)

一方、効率的な数値積分法として、二重指数関数型(DE)公式[3]が提案されている。DE 積分法
では、(7)式のように y = t(x)により任意区間[a,b]の積分を無限区間[-∞,∞]の積分に変数変換し、求積公式として等幅 h の台形公式を適用する。DE 公式では、積分区間の両端に向かって、被積分関数が二重指数関数的に減衰するような変数変換を行う。
dt  x  f  y  dy   f  t  x   dx dx  h  t   x  f t  x 
b



a



i 

i

i

(7)

半無限区間[0,∞]の積分に対する変換関数として３種類の変換が提案されているが、本研究ではマッピングパラメータを導入した(8)-(10)式の形式で適用した。

t  x   exp  sinh x (8)t  x   exp  x  exp   x  (9)t  x   ln  exp  sinh x   1(10)(7)式の求積公式と(8)-(10)式の各変数変換を組合せると、(11)-(13)式の各動径グリッド(DE1,DE2, DE3 グリッド)が得られる。
 F  r  r dr  h   exp 3 sinh x  cosh x F  r 

(11)

2

0

i

i 





0

i 

i

i

2
 F  r  r dr  h  exp 3 xi  3exp  xi    exp   xi  F  ri 





0

i 

2
2
 F  r  r dr  h  ln  exp  sinh xi   1

exp  sinh xi  cosh xi
exp  sinh xi   1(12)

F  ri (13)【計算】正確な積分結果が全電子数となることが分っているため、被積分関数として電子密度を用いた。希ガスを除く第１周期から第４周期までの原子から成る中性の２原子分子種に対して、
TA, MK, DE1, DE2, DE3 の各動径グリッドを適用した数値積分を実行し、積分結果の収束性と精度について性能比較を行った。本研究では、等核２原子分子 A2 及び水素化物 AH と酸化物 AO (A =H, Li-F, Na-Cl, K-Br)を検討した。
基底関数は、H, Li-F, Na-Cl, K, Ca, Ga-Br に対して Pople の 6-31G**基底を適用し、Sc-Zn に対して Ahlrichs の SVP 基底を適用した。分子軌道は、Gaussian09 の M06-DFT 計算による構造最適化の結果を用いた。分子積分は、Becke のファジーセル法により原子積分の和へ分解した。(11)-(13)式の DE 公式における無限和は、適当な範囲に分布したグリッド点に対応する有限和で近似した。
角度積分に対する数値積分近似に起因する誤差を最小にするため、角度グリッドはグリッド点の数が 1202 の Lebedev グリッドを採用した。
動径グリッドのグリッド点の数は 30 から 200 まで 10 ずつ変化させ、積分結果の収束性を確認した。積分精度は、Accuracy = – log|(Approx/Exact) –1|の指標を用いて検証した。ここで、Approx
は数値積分の値を表しており、Exact は正確な値(全電子数)を表している。なお、積分精度の上限(machine precision)は、倍精度の Fortran90 を用いて Accuracy = 15.7 である。
【結果】グリッド点の数が n = 50, 100, 150, 200 の動径グリッドを用いた場合について、A2, AH, AOの各分子種に対する Accuracy の平均値を表に示す。n = 50 では TA グリッド及び MK グリッドがDE グリッドよりも精度のよい結果を与えるが、正確な値への収束は DE グリッドが TA グリッド及び MK グリッドよりも優れている。TA グリッドと MK グリッドで平均 Accuracy の収束が遅いのは、金属化合物に対する収束が非常に悪いためである。結果の詳細は、発表当日に報告する。
表．各動径グリッドによる２原子分子種の電子密度の数値積分に対する積分精度の平均値
n
50
100
150
200

TA

MK

DE1

DE2

DE3

A2

AH

AO

A2

AH

AO

A2

AH

AO

A2

AH

AO

A2

AH

AO

8.0
13.0
15.1
15.4

7.8
11.9
13.2
13.8

7.8
12.5
14.0
14.6

8.7
12.8
14.3
14.8

8.6
12.5
14.1
14.5

8.6
12.6
14.1
14.7

7.3
12.6
15.2
15.5

7.7
13.4
15.3
15.4

7.7
13.5
15.4
15.4

6.7
13.0
15.4
15.3

6.7
13.1
15.4
15.3

6.6
12.9
15.4
15.4

6.6
12.0
15.2
15.4

6.7
12.0
15.2
15.3

6.8
12.0
15.3
15.4

[1] O. Treutler, R. Ahlrichs, J. Chem. Phys., 102, 346-354 (1995).
[2] M. E. Mura, P. J. Knowles, J. Chem. Phys., 104, 9848-9858 (1996).
[3] H. Takahashi, M. Mori, Publ. RIMS, Kyoto Univ., 9, 721-741 (1974).

1P-108
ゼロ磁場分裂定数 D の分子軌道法による計算(阪大院理)

	
  吉村 翔平, 中島 諒, 伊藤 章, 木下 啓二,
川上 貴資, 北河 康隆, 山中 秀介, 奥村 光隆

Ab initio calculations of zero-field splitting parameter D
by the molecular orbital methods(Osaka Univ.)

S.Yoshimura, R.Nakajima, A.Ito, K.Kinosita
T.Kawakami, Y.Kitagawa, S.Yamanaka, M.Okumura【序】	  三重項以上の多重度をもつ分子において観測されるゼロ磁場分裂は、その電子構造を大きく反映す
る重要な物理量であり、分子磁性の分野において実験、理論計算の両面で積極的に研究が行われている。しかし、今日に至るまで、その相互作用パラメーターで異方的なスピン間相互作用に関連づけられる D 値、E 値の定量的な計算手法は確立されておらず、量子化学計算によるアプローチが期待されている。ゼロ磁場分裂は Spin-Spin Coupling(SS)項と Spin-Orbit Coupling(SOC)項に起因し、これらの項は 2 次の摂動論により記述できる。主に金属錯体では SOC 項が、有機分子では SS 項が主要な項であることが多いが、両者の効果がともに無視できない場合もある。
	  このゼロ磁場分裂定数の計算を行ううえで、我々は単分子磁石に注目した。単分子磁石は分子内の磁気異方性により高スピン状態が安定となる物質群であり、分子メモリなどのナノデバイスへの応用が期待されている。この単分子磁石においては、高スピン状態が安定となるためにはゼロ磁場分裂定
数 D が大きな負の値をもつことが実用上必要である。しかし現状では、この D 値の符号を決定する要素に関する知見はあまり得られておらず、系統的な単分子磁石合成の指針は与えられていない。そこ
で、われわれは D 値の符号に関する規則性を見出だすことを念頭に、量子化学計算によるアプローチを行った。
	  今回は、MnII(terpy)X2 錯体についての計算を行った。この系は配位するアニオン(X—)により D 値の符号が変化することが報告されている興味深い系である。

【理論】
	  ゼロ磁場分裂は以下のハミルトニアンで記述される。

D 値への寄与は、SS 項については以下のような McWeeny-Mizuno の式にスピン密度を代入して得られる。

1P-108
SOC 項は以下のような式を用いて得られる。

【計算】	  今回は MnII(terpy)Cl2、MnII(terpy)Br2、 MnII(terpy)(SCN)2 について、D 値の符号に対する分子構
造、及び配位子の違いによる効果を検証するために、種々の計算を行った。計算は密度汎関数法により行い、汎関数は UB3LYP 汎関数を用いた。基底関数は中心金属である Mn(II)については wachters+f、配位子については 6-31++G**、6-311++G**を用いた。また、D 値の計算は、SS 項については得られた UNO のスピン密度を McWeeny-Mizuno の式に代入して求めた。SOC 項については、Neese らの
Coupled-Purturbed(CP)法、及び Pederson、Khanna らの Pederson-Kahnna(PK)法に基づいて行った。全ての量子化学計算には Gaussian09、ORCA2.9.1 を用いた。結果の詳細は当日報告する。

	
 
	
  	
  	
  	
  	
 
	
  	
  	
  	
  	
  	
  	
  	
  	
  	
  	
 

	
  	
 

図 1. Mn(terpy)Br2【参考文献】[1] F.Neese, E.I.Solomon, Inorg. Chem. , 37, 6568, 1998
[2] F.Neese, J.Chem. Phys., 127, 164112, 2007
[3] C.Mantel, C.Baffert, I.Romero, A.Deronzier, J.Pécaut, M.N.Collomb, C.Duboc, Inorg. Chem. 43, 6455,
2004

1P109
HAunCl モデル分子からの Cl 脱離機構に対する理論的研究（阪大院・理）多田 幸平，林 祥生，

坂田 晃平，

北河 康隆， 川上貴資，山中秀介， 奥村 光隆【序】金を粒径 5nm 以下で金属残化物に担持したものは金担持触媒と呼ばれる。金担持触媒には、
例えば-70℃の低温でも一酸化炭素の酸化反応に対して活性を示すといった、他の貴金属担持触媒には見られない特異な触媒活性が見出されている。金を粒径 5nm 以下で金属酸化物に担持するには、中和や水素ガスによる還元により塩素を除去するか前駆体として塩素を含まない金化合物を用いなくてはならない。このような、金担持触媒調製時における塩素による金粒子の凝集促進はよく知られたことであるが、その詳細な機構についての明快な説明は未だなされていない。そして、塩素による金微粒子凝集機構を正しく理解するためには、塩素が除去されることにより凝集がおこらなくなること、塩素が脱離していく機構も調べる必要がある。そのための準備として、
HAunCl モデル分子からの HCl 分子脱離反応を用いて、密度汎関数理論（DFT）に基づく量子化
学計算において度々問題になる、「基底関数依存性」と「汎関数依存性」を検証した。
【計算手法】計算プログラムには Gaussian09 を用いた。金の基底関数依存性を調べるために、
SDD、LANL2DZ、CEP-31G、def2-QZVPD、aug-cc-pVQZ を用い、各々対応した擬ポテンシャ
ルも使用した。塩素の基底関数依存性を調べるためには、6-31+G**、6-31G**、LANL2DZ(d,p)
を用い、LANL2DZ(d,p)に対してのみ、対応する擬ポテンシャルを用いた。交換・相関汎関数依
存性については、PBEPBE、BLYP、PBE0、B3LYP、LC-wPBE、CAM-B3LYP を用いて調べた。
【結果と考察】n=2 の場合のモデル反応を図 1 に示す。

IM1
図1

TS

IM2

HAu2Cl モデル分子から HCl 分子が脱離するモデル反応系。黄が金原子を、緑が塩素原子を、白が水素原子を表している。
まず、金の基底関数依存性について述べる。いずれの基底関数においても図 1 に示した中間体構造（IM1、IM2）及び遷移状態構造（TS）が収束した。得られた活性化エネルギーは表 1 に示した。基底関数依存性は大きくないが、LANL2DZ や CEP-31G よりも SDD の方が def2-QZVPDや aug-cc-pVQZ に近い値を与えることが分かった。
次に、塩素の基底関数依存性について述べる。この場合も同様に、図 1 に示した全ての構造が収束し、構造に対する基底関数依存性は見られなかった。また、表 1 に示した活性化エネルギーについても依存性は見られなかった。

最後に、汎関数依存性について述べる。いずれの交換・相関汎関数を用いても図 1 に示した構造は収束し、構造への汎関数依存性はあらわれなかった。各交換・相関汎関数で得られた構造に対して CCSD(T)/def2-QZVPD レベルで一点計算を行い、それにより見積もられる活性化エネル
ギーEa(CCSD(T))と、DFT により見積もられた Ea(DFT)を比較した。結果は、表 2 に示した。
pure-DFT よりも hybrid-DFT の方が CCSD(T)に近い値を与える傾向があり、特に PBE0 が
CCSD(T)とよく一致した。他の n の場合については、当日発表する。

表 1 HAu2Cl モデル反応系における活性化エネルギーEa の基底関数依存性。水素の基底関数には、6-31+G**を用い、交換・相関汎関数には PBE0 を用いた。

basis set for Au

basis set for Cl

Eaa /kcalmol-1

SDD

6-31+G**

26.92

LANL2DZ

6-31+G**

25.50

CEP-31G

6-31+G**

27.56

def2-QZVPD

6-31+G**

26.87

aug-cc-pVQZ

6-31+G**

26.73

SDD

6-31+G**

26.92

SDD

6-31G**

26.72

SDD

LANL2DZ(d,p)26.91a) Ea = E(TS) - E(IM1)

表 2 HAu2Cl モデル反応系における活性化エネルギーEa の交換・相関汎関数依存性。水素と塩素の基底には 6-31+G**を用いた。金の基底には DFT の場合に SDD を、CCSD(T)の場合に
def2-QZVPD を使用した。エネルギーの単位は全て kcalmol-1。

交換・相関汎関数

Ea(DFT)

Ea(CCSD(T))

Ea(DFT) - Ea(CCSD(T))

PBEPBE

21.71

27.04

-5.33

BLYP

20.66

26.89

-6.23

PBE0

26.92

26.63

0.29

B3LYP

24.93

26.79

-1.86

LC-wPBE

35.89

26.42

9.47

CAM-B3LYP

29.81

26.51

3.30

1P-110

メニーコアプロセッサによる量子化学とその応用（名大エコトピア1，㈱クロスアビリティ2）安田 耕二1，古川 祐貴2，古賀 良太2

Quantum Chemistry on Many-core Processors（Nagoua Univ.1，X-Ability Co.Ltd.2）Koji Yasuda1，Yuki Furukawa2，Ryota Koga2【序】微細加工技術に伴って計算機は高速化してきたが、消費電力が問題になりつつある。現在
では 1 語をメモリーに転送すると、数百演算分の電力を消費する。つまり電力面からはメモリーを読み書きしないアルゴリズムが有利になる。その結果量子化学計算では、一個の二電子積分はおよそ数十演算で計算できるので、direct SCF が電力面から望ましくなる。更に積分計算アルゴリ
ズムは、計算量よりメモリーの読み書き回数で優劣が決まる。また計算機も、単純なコアを多数集積したメニーコア型が、従来型より一桁電力効率が良いため、主流になると思われる。我々はメニーコア型計算機(GPU)での量子化学計算法を研究してきたが、今回はメモリー読み書きを減らすアルゴリズムと、d 基底関数への実装を報告する。
【方法】漸化式による二電子積分計算では、多数の中間積分ができる。これがレジスタ上に保持できないと、キャッシュやメモリーに待避される。このレジスタ溢れを押える事が、速度向上と電力削減に重要である。メモリー操作回数で漸化式の優劣を測る指標として、MOpコストが知られている。これは漸化式の右辺に現れる異なる中間積分の総数を示す。MOpコストは、レジスタ上に必要な中間積分を全て保存できる場合、レジスタ溢れや電力消費の良い指標になる。が現在のCPUはもっと少ないレジスタしか持たない。そこで我々は、レジスタ溢れを押える積分計算法
について検討した。具体的にはMcMurchie- Davidsonの一中心積分[r](m)から、静電ポテンシャル行
列要素をdirect法で計算する際、①どの順序で積分を計算するか ②どう並列化するか、について検討した。
① 一中心積分[r] (m)をBoys関数[0] (m)に還元するには、最小の非ゼロ指数を持つ軸（x, y, zどれか）方向の漸化式を使えば良い。ある中間積分は、漸化式の右辺に複数回現れる。計算した中間積分
はレジスタから溢れる前に、早く使い切るのが望ましい。そこでどの順序で中間積分を計算すると、中間積分を保持するレジスタ数を最小化できるか調べた。それには中間積分のデータ依存関係を木構造に表し、どの部分木から評価するのが良いか考えればよい。これは、部分的にはコンパイラが最適化時にやるのだが、この最適化はNP困難な問題で解けない。我々はこの最適解を総当りで調べた。
② 二電子積分は普通 shell 四つ組単位で計算される。そこで複数の shell 四つ組を同時に計算して、並列化やベクトル化を行う。この時静電ポテンシャル行列要素をどこに保持するかも問題である。d 関数では 35 個の角運動量成分があるが、レジスタに置くには数が多く、メモリーに置くと読み書きが増える。一中心積分計算は x, y, z 方向に対称的なので、我々は 3 個のスレッドで１
つの shell 四つ組の計算をし、必要なレジスタ数を減らす方法を検討した。また静電ポテンシャル
行列要素も 3 個のスレッドのレジスタに分散して割り付けた。必要な積分を他のスレッドが持つ
場合、再計算した。
【結果】中間積分を保存するのに必要な最小レジスタ数は、中間積分の総数よりかなり小さく、
例えばL=4, 8 でそれそれ 6, 12 だった。適切な計算順序ではレジスタ数は少なくてすむ可能性を示す。次にmが小さい順に[r] (m)を計算したソースコードと、最適順序で計算したソースコードをnvcc
でコンパイルし、使用レジスタ数を調べた。その結果角運動量Lが小さいなら差はなく、Lが大き
いと差があった。このレジスタ数には、静電ポテンシャル行列要素やループカウンターなども含んでいる。また 3 スレッドでshell四つ組を計算する場合、Lが大きくてもレジスタ数はかなり少なくなった。また、これらのソースコードで 1024 個のshell pair間の二電子積分から静電ポテンシャ
ル行列要素をGPUを用いて計算した時の性能を測定した。使用レジスタ数 63 以下では、理論性能の約 5 割を示し、それ以上では約 2 割に性能は減少した。
次にこの積分計算法を Gaussian と GAMESS-FMO に組み込んだ。
Gaussian：Linda で並列化した上で、各 node の計算は

表１：SCF10 反復にかかった時間（秒）OpenMP と GPU で高速化した。次の問題があった。

original

+GPU

21.8

93.6

2440.228.8

二電子積分(GPU)

-

55.4

交換相関項(CPU)再計算した。(2) Gaussian09 の高速多重極計算ルーチンには、 交換相関項(GPU)

211.5

1.2

-

80.8

50.7

27.8

2744.4313.9(1) GPUで各nodeの計算速度を 10 倍程度加速した結果、

高速多重極

Lindaの通信時間が無視できなくなった。Gaussian09 では、 二電子積分(CPU)
二電子積分計算時に、全てのサブルーチン引数を通信で送っている。そこでできるだけ通信をやめ、各nodeで引数を競合と思われるバグがあり、OpenMP並列で大きな計算をするとたまに間違える。(3) GPU計算がまれにunspecified

通信
計 (SCF10 反復）

launch failure（原因不明のでエラー）で終了する事があった。この再現困難なエラーはGPUのスタック領域不足が原因のようだ。(4) Second order SCF(Newton法)をGPUで加速する場合、連立方程式の反復解法が数値誤差の影響を受け易い事が分かった。そこでDIIS法が失敗する場合は、scaled steepest decent法を用いた。水和白金クラスター[Pt13(H2O)15]、基底関数はD95 とStuttgart/Dresden ECPを用い（801 AOs）、2 node×6 threadsで計算した時間を表 1 に示す。全体で 9 倍程度加速された事、Linda並列化されていない高速多重極計算や、通信時間が無視できなくなっている事が分かる。
GAMESS-FMO：フラグメント分子軌道法(FMO 法)の環境静電ポテンシャル(ESP)の計算の加速を行った。モデル分子はヒトインシュリン(PDBID:2HIU)、CPU は Intel Core i7-3930K 3.2GHz (6core)、GPU は NVIDIA GeForce GTX680, コンパイラには Intel FORTRAN/C++ Compiler と CUDA4.0を用いた。この計算には、クロスアビリティが販売している GPU 計算エンジンモジュール XACUDA-QM [1] を用いた。

表 2 ： ESP 部 分 と 計 算 全 体 の 加 速 率
表 2 に種々の条件での計算全体とESP部分のみの加 （GAMESS を 1 とする）

速率をまとめた。どの場合も全体でGAMESSを 2.6
-6

ESP

全体

～4.4 倍加速できた。また、エネルギー値の差は 10

6-31G, FMO2, 一点計算

14

4.1

a.u.以下であり、十分小さかった。また同様の計算

6-31G*, FMO2, 一点計算

17

4.4

をGeForce GTX580 GPUで行ったが、加速率はほぼ

6-31G, FMO3, 一点計算

20

3.7

同じだった。GTX680 の理論性能がGTX580 の 2 倍

6-31G, FMO2, 勾配計算

11

2.6

である事を考えると、改善の余地がまだある事が分かる。詳しい解析は当日報告する。
[1] Furukawa, Y., Koga, R., Yasuda, K. JSST 2011 Tokyo. (2011)

1P111Elongation 法によるポルフィリンワイヤーの効率的物性計算(九大院・総理工 1，SCNU 2，JST-CREST 3)
阿部 真也 1，Liu Kai 1，折本 裕一 1，Gu Feng Long 2,3，青木 百合子 1,3Theoretical study on metalloporphyrin array by elongation method(Kyushu Univ.1, SCNU2, JST-CREST3)
Shinya Abe1, Liu Kai1, Yuuichi Orimoto1, Gu Feng Long2,3, Yuriko Aoki1,3【序】ポルフィリンの金属錯体を含む骨格は生体系において光合成の促進、および電子伝達に関して重要な役割を果たしている。そのような優れた電子的・光化学的特性を示すことから新規機
能性材料を目指したポルフィリン多量体の合成が盛んに行われており、ポルフィリン同士がメゾ位において直接結合した多量体についてはこれまでに最長で亜鉛(Ⅱ)ポルフィリン 1024 量体が報告されている。1) また、ポルフィリンは大きく広がった共役系π電子と 1 次元に非局在化した電子の存在から、3 次の高分子非線形光学(NLO)材料として注目されている。一方当研究室では、局在化と相互作用を繰り返しながら計算機上で高分子の重合反応を模倣することにより電子状態を効率的に得るための新規な手法-Elongation(ELG)法-を開発してきた。しかし、系全体に非局在化したπ軌道等を有する系については、領域局在化軌道を作る際に攻撃フラグメントに近い領域
へと十分に局在化されないことにより、計算誤差が大きくなることが問題となっていた。この問
題点を踏まえて ELG 法を改良し、様々なポルフィリンアレイに適用することで、非線形光学材料
および導電性材料の設計に向けての物性計算を試みた。
【理論】ELG 法とは、適当な大きさの出発クラスターに対して、任意のフラグメントを付加して系を次々と伸長させながら高分子鎖の電子状態を効率的に求める方法である。具体的には、まず出発クラスターの Hartree-Fock 方程式を解いて、系全体に広がった正準分子軌道(CMO)を得る。
次に、出発クラスターに対して Frozen(A)領域、Active(B)領域を定義し、CMO に対するユニタ
リー変換によって、A 領域、B 領域それぞれに局在化した領域局在化分子軌道(RLMO)をつくる。
続いて、新しいモノマーを B 領域の末端から付加させて、高分子鎖を伸長させる。この際、攻撃分子を B 領域に局在化している RLMO とのみ相互作用させることにより、対角化させる Fock 行
列の次元を減らすことができる。A 領域と攻撃分子の間の距離は十分離れているため相互作用せず、計算結果には影響しない。その結果、高分子鎖を伸長しても対角化に要する次元の大きさは常に一定となり、全系をまともに扱う従来の計算手法と比べて、極めて効率良く電子状態を求めることができる。しかしながら、HOMO、LUMO を構成するπ軌道が分子鎖全体にかけて強く非局在化した系においては特定の領域への局在化が満足に達成されず、本質的に局在化できない
いくつかの軌道が Frozen 領域に残ってしまう。そこで我々は、ELG 法において以下のような
Orbital Shift(OS)法を導入した(OS-ELG 法)。2) ここで、領域局在化において Frozen 領域中の軌道を拾い出して Active 軌道として付加分子と相互作用させるかどうかを選択するための Frozen
RLMO と Active RLMO の重なり積分の値は次のように定義される。

Qij = φiFrozen RLMO ( A) φ jActive RLMO ( B)

設 定 し た あ る 閾 値 よ り も 大 き な Qij を も つ
Frozen RLMO が存在するとき、それらの軌道はActive RLMO に帰属され、付加分子との相互作用に含まれる。この手法を全ての伸長ステップ
を通して続けることで、非局在化軌道は常にActive 空間へとシフトする。さらに各ステップ
において分子(超)分極率を得るために、鎖長方向にかけた電場による摂動を含めた全系の効率的な 計 算 が 可 能 な Elongation-Finite Field(ELG-FF)法 3)をポルフィリンアレイに適用した。
ここで、分子(超)分極率は各電場で得られた全エ
ネルギーを電場(E)で数値微分し、E =0 における微分係数から求められる。(Fig.1)【結果と考察】ポルフィリンの隣り合う平面ユニット間の二面角を 45°とした系に対し、ELG法により 1 次元方向に伸長させて RHF/6-31G**
レ ベ ル で 計 算 し た 場 合 の 、 旧 方 法 (ELG) と
OS-ELG 法(Qij =10-4)それぞれの全エネルギー

Fig. 1Elongation-FF method

誤差の比較を Table.1 に示す。OS-ELG 法では 3~4 桁ほど精度が向上しており、従来法との誤差は 10-9 a.u./atom と極めて小さい。OS-ELG 法によって大きなπ電子系でも非局在化したπ軌道に由来する誤差が大幅に改善されることが分かった。ELG-FF 法から得られた鎖長方向の 3 次の超分極率 γzzzz とその各伸長ステップ毎の増分 Δγzzzz を Fig.2 に示す。従来法(CONV)に対しては全エネルギーの高次微分に起因する数値的不安定性が見られるが、OS-ELG に対しては SCF 計算の収束性も良好で、なおかつ γ 値においても旧 ELG 法と比較してより全系計算と一致する結果が得られている。また、非局在化したπ軌道の NLO 特性への影響の大きさを示しているとも言える。
今後は平面縮環型ポルフィリンアレイや中心金属を含む系についても本手法の有効性を検討し、バンド構造から導電性と分子構造の関係についても調べていく。

Table. 1 全エネルギーの誤差

Fig. 2γzzzz、Δγzzzz の鎖長依存性【参考文献】1) A. Osuka, et.al. Chem.Eur.J.,11,3389-3404 (2005)
2) Y. Aoki, O. Loboda, K. Liu, M. A. Makowski, and F.L. Gu, Theor. Chem. Acc., 130, 595-608 (2011)
3) F. L. Gu, Y. Aoki, A. Imamura, D. M. Bishop, and B. Kirtman, Mol. Phys., 101, 1487-1494 (2003)

1P112

Rigged QED 数値シミュレーションによる電子スピンダイナミクスの理論的研究（京大院工）福田 将大, 市川 和秀, 立花 明知

Theoretical study of electron spin dynamics
by Rigged QED simulation（Kyoto Univ.）Masahiro Fukuda, Kazuhide Ichikawa, Akitomo Tachibana
fukuda.masahiro.57w@st.kyoto-u.ac.jp

原子・分子系のような微小領域において、光子と電子のスピンダイナミクスの数値シ
ミュレーションを行うには、場の量子論に基づいた計算を行う必要がある。本研究では、場の量子論的に原子と光子の相互作用を扱うことが可能である Rigged QED (Quantum

ElectroDynamics) [1] に基づく数値シミュレーションにより、時々刻々と変化する電子のスピン角運動量密度を可視化し、そのダイナミクスの描像から旋光性について議論する。
電子のスピン角運動量密度の運動方程式は次のように表わされる [2]。

∂ˆˆ⃗se (x) = ⃗tˆe (x) + ζ⃗e (x).
∂t(1)

式 (1) の右辺第一項はスピントルク密度であり [2]、t̂ke (x) = −εlnk τ̂eΠ ln (x)

(2)

Π ln
Π ln
= −εlnk τ̂eN(x) − εlnk τ̂eA
(x)

(3)

Π ln
Π ln
(x) + τ̂eA(x) はストレス
ここで、εlnk はレビ・チビタテンソルであり、τ̂eΠ ln (x) = τ̂eNテンソルを表す。

)
i~c ( †
ψ̂ (x)γ 0 γ l ∂k ψ̂(x) − (∂k ψ̂ † (x))γ 0 γ l ψ̂(x)2)
Ze e ( †Π klτ̂eA(x) = −
ψ̂ (x)γ 0 γ l Âk (x)ψ̂(x) + ψ̂ † (x)Âk (x)γ 0 γ l ψ̂(x)
2
Π kl
(x) =
τ̂eN

(4)
(5)

⃗ˆ
ψ̂(x) は場の演算子、A(x)はベクトルポテンシャル、γ はガンマ行列を表す。
式 (1) の右辺第二項はツェータ力密度である [2]。
(
)
k1
k
k
ˆ
ζ̂e (x) = −c∂k ψ̄ (x)γ ~Σ ψ̂(x) ;2(no sum over k)(6)⃗ はパウリ行列を対角成分に持つ 4×4 行列である。
Σ⃗ˆ⃗ˆA (x) + A⃗ˆrad (x) という遅延ポテンシャル A⃗ˆA (x) と外ベクトルポテンシャルは A(x)=A⃗ˆrad (x) の和で表わされる。遅延ポテンシャルが無視できるほど小さいとする部光子場 Aと、スピントルク密度は次のように表わされる。
Π ln
Π ln
t̂ke (x) ≈ −εlnk τ̂eN(x) − εlnk τ̂eArad
(x)

(7)

ただし、

1 l
ĵ (x)Âk (x)c e [ rad]ĵek (x) = Ze e c ψ̂ † (x)γ 0 γ k ψ̂(x)Π kl
τ̂eArad
(x) =

(8)
(9)

式 (7) の右辺第二項は、電子と光子の相互作用によって光子から電子スピンに与えられるトルクを表している。
式 (1)(7) から次のような理論的な説明がなされる。外部光子場の無い定常状態の電子
スピンは局所的にスピントルクとツェータ力が拮抗することにより安定な状態を保っている [3]。しかし、電子が外部からの光子と相互作用する場合には光子から受けるトルク
によりそのバランスは崩れ、電子のスピン角運動量密度は時々刻々と変化する。一方で、光子が電子スピンに与えるトルクの反作用として光子にもトルクが働く。これが旋光性という性質の起源であると考えることができる。(もっとも、場の演算子 ψ̂(x) で表わされる電子状態によって決まる選択則から遷移が禁止される場合には旋光性は生じない。）
以上のような原子・分子系で起こる局所的な電子と光子のスピンダイナミクスを Rigged

QED 数値シミュレーション [4] によって調べ、その描像を発表で示す予定である。

参考文献[1] A. Tachibana, in Fundamental World of Quantum Chemistry, A Tribute to theMemory of Per-Olov Löwdin, ed. by E. J. Brändas and E. S. Kryachko, (Kluwer
Academic, Dordrecht, 2003), Vol. 2, p. 211.
[2] A. Tachibana, J. Mol. Struct. (THEOCHEM), 943, 138 (2010).
[3] T. Hara, M. Senami, A. Tachibana, Phys. Lett. A 376, 1434-1441 (2012).
[4] K. Ichikawa, M. Fukuda and A. Tachibana, Int. J. Quantum Chem., published
online (2012). [ DOI: 10.1002/qua.24087 ]1P113
反応経路探索法 GRRM における制限探索オプション LADD の効能の評価(和歌山大学 1,豊田理化学研究所 2) 高田谷吉智 1, 山門英雄 1,大野公一 2
Parformance of the GRRM-large ADD following option for explorationof chemical reaction pathways(Wakayama Univ.1,Toyota Phys. Chem. Res. Inst.2)
Yoshitomo Kodaya1,Hideo Yamakado1,Koichi Ohno2【序】 従来 4 原子以上では不可能とされてきた化学反応経路の全面自動探索は、超球面探索法[1]によって可能になったが、原子数が増えると探索時間が非常に長くなる。非調和下方歪み(ADD)の大きい経路を優先して探索する Large ADD following (lADDF)法を用いると限定的な探索にはなるが、化学的に重要な低エネルギー領域が優先的かつ効率的に探索される。そこで、lADDF 法のオプションを変化させ、その効能を解析した。
【方法】 GRRM 11 プログラム[2] を用い、オプシ

100ョン変数の NRUN や LADD を変化させ、効能を評価

BCNOS分子した。また、量子化学計算には Gaussian プログラム
路を何番目まで辿るかの指定であり、ADD は、振動固有値の平方根でスケールした基準座標で記述した
調和ポテンシャルと実際のポテンシャルとの差であ△:LADD=5
Number of EQs（g03, g09）を使用した。LADD は ADD の大きい経50□:LADD=3○:LADD=1る。また、NRUN は、乱数を用い原子をある範囲内
でばらまいて構造最適化をするときに、自動発生させる初期構造の数の指定である。全面探索結果が既知の0

BCNOS 分子[3]について、LADD と NRUN を変化さ

0

20

40

NRUN

せ、探索される平衡構造(EQ)と遷移構造(TS)の数、

図 1. LADD=1, 3, 5 のときの各 NRUN

LADD や NRUN を変化させたときの EQ のエネルギに対する探索された EQ の数。

ー分布、全面探索と LADD の制限を付けたときの計

200

算時間の比率を調べた。

BCNOS分子5, 10, 20, 30, 50 としたときの EQ と TS の数をそれぞれ図 1,2 に示した。LADD=1 についてはばらつきを減らすため 10 回分の平均を示した。NRUN の数を増やすと見つかる EQ、TS の数は増えるが LADD=3, 5 で

Number of TSs【結果・考察】 LADD=1, 3, 5 について NRUN= 1, 2,△:LADD=5100□:LADD=3

は NRUN=30 から 50 に増やしても見つかる数にほぼ変化が見られず、飽和していく傾向がみられる。
図 3, 4 は LADD=3, 5 について各 NRUN に対して
EQ のエネルギー分布を示したものである。NRUN=2
以降は増加分だけを示し、最安定構造である○:LADD=1
0
0

20

40

NRUN
図 2. LADD=1, 3, 5 のときの各 NRUNに対する探索された TS の数。

BCNOS分子

BCNOS分子500

LADD=3△:NRUN=1▲:NRUN=2～50
(増加分)

0
0

20

Energy / kJ・mol

Energy / kJ・mol

−1

−1

1000

1000

500

LADD=5△:NRUN=1▲:NRUN=2～50
(増加分)

0

40

0

20

NRUN

図 3. LADD=3 についての NRUN に対

40

NRUN

図 4. LADD=5 についての NRUN に対するエネルギー分布図するエネルギー分布図

SCNBO(線形分子)を 0 kJ/mol とする相対エネルギーBCNOS分子

LADD=3 では 300 kJ/mol 以下の安定な構造の大部分が、
LADD=5 では 500 kJ/mol 以下の構造はすべてが、探索されている。NRUN を増やすと、EQ は高エネル
ギー領域に新たに追加される傾向が見られた。

Energy / kJ・mol

だけで効率的に探索できている。NRUN=1 で、−1で表示した。低エネルギー領域については、NRUN=1
1000

500

NRUN=1△:LADD=1▲:LADD=2～10(増加分)

LADD=5 では、NRUN= 30 以降で新たに見つかった
EQ は僅かであった。

0

図 5 は NRUN を一定にし、LADD を 1 から 10 まで変化させたエネルギー分布図である。LADD を増加させると、新たに見つかった EQ は徐々に高エネルギー側にシフトしていく傾向にある。

0

5

10

LADD

図 5.

NRUN=1 について各 LADD に対するエネルギー分布図

LADD=8 付近からは、一番高いエネルギーの EQ のエネルギーにほぼ変化がなく、飽和していく傾向が見られた。図 5 のエネルギー分布の下端から、指定された高さのエネルギーまで効率的に調べるために必要な LADD の大きさを知ることができる。たとえば、300 kJ/mol 以下の低エネルギー領域については、LADD=3 でもすべての EQ が探索されることがわかる。
【結論】 NRUN を増やせば探索される EQ や TS の個数は増えるが、NRUN=30 以降はあまり増加が認められない。調べたい低エネルギー領域の上限を指定すれば、それ以下の EQ をすべて探索するのに必要な LADD の大きさを、推定できることがわかった。低エネルギー領域を優先す
る探索については、NRUN は 1 でも十分に機能することが認められた。
[1] K. Ohno and S. Maeda, Chem. Phys. Lett. 348, 277 (2004); S. Maeda and K. Ohno, J. Phys.
Chem. A109, 5724 (2005); K. Ohno and S. Meda, J. Phys. Chem. A110, 8933 (2006).
[2] 大野公一, 長田有人, 前田理, 諸熊奎治, 第 14 回理論化学討論会, 岡山 (2011), 2D1b.
[3] Koichi Ohno and Yuto Osada, "Advances in the Theory of Quantum Systems in Chemistry andPhysics", Springer, 381-394 (2012).

1P114
長距離間相互作用補正密度汎関数法による分子の動的超分極率算定（豊橋技術科学大学）石丸広騎、浜田信次、関野 秀男[An evaluation of dynamic hyperpolarizability by Long range
corrected Density Functional Theory method]（Toyohashi University of Technology）Hiroki Ishimaru,Shinji Hmada,Hideo Sekino【研究背景・目的】
物質の光学物性を調べることは重要であり，光学物性を知ることは工業的利用を目的とした光機能材料の開発に大いに役立つ．物質の光物性を調べるためには分子の分極率および超分極率を解析することが必要であるが，特に物質の高次の応答である超分極率の定量測定は困難であり，計算量子化学の貢献が期待される．
Ref１では小さな分子において実測値と近い値が得られており，定量算定には電子相関の考慮が重要であることが明らかだが，それらを考慮した方法として， Coupled-cluster singles anddoubles（CCSD）法や Density Functional Theory（DFT）法が挙げられる．CCSD 法は厳密な手法であり小さな分子での実測を再現できる．しかし，CCSD 法では計算コストが非常に高く，分子のサイズが大きくなる場合には膨大な時間がかかってしまう．そこで，本研究では HF 近似と同程度のコストで電子相関を考慮している，DFT 法を採用した．しかし，DFT 法での大きな分子の分極率
や超分極率の算定には問題があることが知られているため，我々は長距離間相互作用補正法を使ってその問題に対応した．また，超分極率の実測はレーザー入力による定周波数条件下で行われるため，超分極率の周波数依存性を算定する．

【分極率および超分極率】
原子や分子の電子雲が外部電場によって歪められ，歪められた電荷分布と元の状態との偏差を分極率という．
静電場𝜀0 から外部電場𝜀 ±𝑖𝜔𝑡 の影響を受けると分子分極は，
1
1μ(ωt) = 𝜇0 + α𝑒 ±𝑖𝜔𝑡 + 𝛽𝜀0 𝑒 ±𝑖𝜔𝑡 + 𝛾𝜀02 𝑒 ±𝑖𝜔𝑡 + ⋯2!
3!
と展開される．ここで，μを双極子モーメント，𝜇0 を永久双極子，αを分極率，βおよびγを第１および第２超分極率といい，電子の非線形応答に対応し強いコヒーレントな外場でのみ実現する．
本研究では，動的電場に対する波動関数の線形応答を求め，それにより，分極率を解析的に算出した．また，超分極率は Finite-Field（FF）法により更に静電場をかけ，分極率の数値微分により求めた．

【密度汎関数法】
密度汎関数法とは密度汎関数理論に基づく電子状態計算法のことであり，DFT 法は全電子エネ
ルギーを電子密度の汎関数として記述する．HF 法の交換積分項を密度汎関数理論に基づいて電子密度の汎関数に置き換えることで得られる Kohn-Sham 方程式による解法である KS-DFT 法が，現在の DFT 法の主流となっている．サイズの大きい分子の応答を求めるような場合，局所密度近似(LDA)による汎関数を用いた KS-DFT 法の挙動が，破綻をきたす問題が多々報告されている．そのため最新の理論手法では，交換ポテンシャルを距離に応じた割合で取り込み，汎用の交換相関汎関数に補正を行う長距離間相互作用補正[２]を DFT 法に適応した手法が注目されている．

【長距離間相互作用補正】
DFT 法における交換汎関数の長距離での振る舞いを改善するため，長距離間相互作用補正（Long-Range Correction:LC）法が提案された．以下に式を示す．

1
1 − erf⁡(𝜇𝑟12 ) erf⁡(𝜇𝑟12 )
=
+
𝑟12
𝑟12
𝑟12
この補正では，２電子間距離を短距離成分と長距離成分に分割する．𝑟12 は２電子間の距離，μは補正の割合を決める係数である．これにより交換相互作用を算定するのだが第１項は DFT 法により計算され，第２項は HF 法により計算する．この補正により，短距離成分は DFT 法により計算され，２電子間の距離が大きくなるほど，HF 法により計算される割合が大きくなることで，長距離の相互作用が改善される．

【検証方法・結果】
分 子 分 極 率 計 算 で は Coupled Perturbed Hartree-FocK(CPHF) 法 や Local Density
Approximation(LDA)，LC を汎関数 BLYP に適応した LC-BLYP 法を用いた Coupled Perturbed
Kohn-Sham（CPKF）法により，動的な分極率を算定した．第 2 超分極率の算定には前述の FF 法を用いた．以下に１，３，５―ヘキサトリエン（𝐶6𝐻８ ）の結果を示す．基底関数は cc-pVDZ を用いて，計算プログラムパッケージは，NWChem を使用した．
𝐶6𝐻８ や𝐶20 𝐻22 の動的な分極率および超分極率をそれぞれ図１，図２に示す．

図 1．𝐶6 𝐻８ の動的な分極率

図 2．𝐶6 𝐻８ の動的な第２超分極率【考察】
図 1 の分極率について，汎関数として LDA を使用した DFT と HF 近似を比較すると，周波数の低いところで分極率が発散していることがわかる．一方で LC 法を適応した LC-BLYP は LDA よりも
HF 近似に近づいており，過大評価が抑制されていることがわかる．図２の第２超分極率についても同様であり，LDA による過大評価を LC-BLYP は抑制していることがわかる．また，第 2 超分極率では DFT の破綻がより激しく現れることも観測される．より長い分子についの超分極率算定の数値結果は発表時に報告する．

【まとめ】
以上の考察より，LC 法を適応したことで DFT 法における交換汎関数の長距離における相互作用の不適切な振る舞いが改善されたことが分かった．また，第 2 超分極率の周波数依存性をはじめて算定した．
今後は CCSD やその他の汎関数でも計算し，LC-BLYP の有用性を考慮したうえで分子のサイズが大きい場合の分極率および超分極率について解析していきたい．

【参考文献】[１]Hideo Sekino and Rodeny J.Bartlett，J.Chem.Phys.94(1991)3665， Hyperpolarizabilities
of molecules with frequency dependence and electoron correlation[２]松本啓紀，関野秀男，量子科学的手法を用いた高分子光学物性量算定に関する研究（豊橋技術科学大学修士論文(2007)）

1P115有限長 pentaheptite ナノチューブの量子化学的研究（明治薬大）溝口 則幸

Quantum Chemical Study of Finite-length Pentaheptite Nanotubes(Meiji Pharmaceutical University) Noriyuki Mizoguchi
1 序論

有限長のアームチェア炭素ナノチューブにおける HOMO-LUMO ギャップは、6 員環の層の数（ナノチューブの長さ）の増大と共に周期 3 で振動することが知られている。我々は、先に、符号を考慮したケクレ構造の総数(ASC)を使って、アームチェア炭素ナノチューブは 6 員環の層の数の増大と共に周期 3 で非結合分子軌道を持つことを示し、これが HOMO-LUMO ギャップの振動が生じる原因であるあることを示した。[1] 炭素ナノチューブのすべての 6 員環を 5 員環と７員環で置き換えたナノチューブは pentaheptite ナノチューブといわれる。π電子の 6 員環のネットワー
クから、５員環と７員環のネットワークへの変換が、有限長のナノチューブの HOMO-LUMO ギ
ャップの振動にどのような影響を与えるかは興味ある問題である。ここでは、ナノチューブの長さを伸ばしたときに、pentaheptite ナノチューブの電子状態がどうのように変化するかを調べた。
計算には半経験的方法の一つである PM3 を用いて、制限無しの構造最適化を行った。

2(7,7)ｍ pentaheptite ナノチューブの HOMO、LUMO、 HOMO-LUMO ギャップ
図 1 に(7,7)２アームチェア炭素ナノチューブと(7,7)2 と(7,7)3 pentaheptite ナノチューブを示した。

ここで、(7,7)ｍ の添え字ｍは六員環の層の数（ナノチューブの長さ）を表す。pentaheptite ナノチ
ューブの端はアームチェア型とジグザグ型の２種類がある。図 1 に示した pentaheptite ナノチュー
ブの一つの端は常にアームチェアであり、もう一方の端は層の数ｍが偶数のときは、アームチェアとなり、層の数ｍが奇数のときは、ジグザグとなる。

図 1 (7,7)２アームチェア炭素ナノチューブと(7,7)２と(7,7)3 pentaheptite ナノチューブ
5 員環と 7 員環の導入はナノチューブを不安定化させると予想される。実際、筒の長さにかかわらず、pentaheptite ナノチューブのエネルギーの計算値は対応するナノチューブのエネルギーの計算値よりも高い。つまり、pentaheptite ナノチューブは対応するナノチューブよりもエネルギー的に不安定であることがわかった。
図 2 に(7,7)m pentaheptite ナノチューブの HOMO と LUMO ののエネルギーのチューブ長依存性を示した。(7,7)m pentaheptite ナノチューブの HOMO のエネルギーは長さの増大と共に、周期２で振動しながら、増大する。 (7,7)m pentaheptite ナノチューブの LUMO のエネルギーは周期２で振動しながら、減少している。HOMO のエネルギーは層の数ｍが偶数のときは、より低い値をとり、層の数ｍが奇数のときは、より高い値をとる。ジグザグ型エッジを持つグラフェンやナノチュー
ブの HOMO のエネルギーはより高くなることが知られている（エッジにより効果）。上にみたHOMO のエネルギーの周期２の振動は、pentaheptite ナノチューブにおけるエッジ効果によるもの
であると説明できる。一方、LUMO のエネルギーは層の数ｍが偶数のときは、より低い値をとり、層の数ｍが奇数のときは、より高い値をとる。

図2(7,7)m pentaheptite ナノチューブの HOMO と LUMO のエネルギーのチューブ長依存性

図 3 に(7,7)m pentaheptite ナノチューブの HOMO-LUMO ギャップのチューブ長依存性を示した。

図 3 (7,7)m pentaheptite ナノチューブの HOMO-LUMO ギャップのチューブ長依存性(7,7)m pentaheptite 型ナノチューブの HOMO-LUMO ギャップはチューブの層の数（長さ）の増加
と共に、非常に小さな振幅であるが周期２で振動しながら、減少している。この振幅の減少は奇数員環の導入により、分子軌道に対する対定理が成り立たないことによるものである。
文献
3p-8

1 溝口 則幸,第 40 回記念フラーレン・ナノチューブ総合シンポジウム, 2011/３,名古屋,

1P116

量子ウォークによる時間依存シュレディンガー方程式の計算(豊橋技術科学大学大学院 1 , 理科学研究所 2)
川畑雅之 1 , 関野秀男 1 , 浜田信次 2Abstract(Toyohashi University of Technology1 , RIKEN2)
Masayuki Kawahata1, Hideo Sekino1 , Shinji Hamada2【序】量子力学の分野において、超微視的粒子の状態は波動関数と呼ばれる概念を用いて記述され、波動関数は量子力学の基本方程式である時間依存シュレディンガー方程式(Time Dependent Schrödinger Equation : TDSE)に従って時間発展する。ポテンシ
ャルを考慮した計算や多体問題など、非常に難解で膨大な時間を要する計算量子力学
において、基本方程式である TDSE をより効率的に解くことは量子力学分野全体において大きな利益となり得る。本研究では、TDSE 計算の高精度化・高速化を目的とし、TDSE の新たな解法である量子ウォーク(Quantum Walk : QW)による計算をシ
ミュレーションし、解析する。
＜時間依存シュレディンガー方程式 : TDSE＞

i
【原理】


 (t )  Hˆ  (t )t量子ウォークは古典ランダムウォークの量子版といわれ、2000 年頃より量子コン
ピュータ関連の話題として本格的に研究され始めた。極限分布などのいくつかの定理が導出され、次第にその性質が明らかになってきている。
[QW の性質]
 量子ウォークを粒子の運動だと考えるとき、その粒子は量子ビットを持ち、量子ビットは複数の複素数を内包する。
 量子ウォークのダイナミクスはユニタリー行列 U で定義される。
 カイラリティが存在し、行列 U をカイラリティで分けた行列を定義する。

図 1-1.古典ランダムウォーク
図 1-2.量子ウォーク
図 1.古典ランダムウォークと量子ランダムウォークのイメージ

【実験】QW では、量子ビットに作用させるユニタリー行列を変化させることによってさま
ざまな発展を得ることが出来る。本研究では、QW 計算を行うプログラムを作成し、ユニタリー行列を変化させて得られる２つの QW について計算を行った。最も有名な
QW であるアダマールウォーク(Hadamard Walk : HW)と、TDSE の解を与えるシュ
レディンガーウォーク(Schrödinger Walk : SW)の２つである。初期分布関数はデルタ関数またはガウシアンとした。
また、従来 TDSE の一般的な解法とされてきた差分法によるシミュレーションを行うプログラムも作成し、その計算結果を SW の計算結果と比較した。
すべての計算は一次元と二次元の格子で行った。
【結果と考察】初期分布をデルタ関数とした場合、QW の計算結果はすべて振動した。それに対してガウシアンを初期分布とした場合、ガウシアンの形を保ちつつ、HW は移流方程式、
SW は TDSE に相当する結果を与えた。
また、差分法による TDSE と QW による TDSE の比較では、1 次元 TDSE においての計算結果は高精度で一致した。2 次元 TDSE おいては、全く別の発展を行う結果が出たが、これは量子ウォークを１次元から２次元へと拡張する過程において、時間
発展の振る舞いを左右するユニタリー行列の設定に任意性があるためだと考えられる。

図 2.１次元 HW(delta)

図 3.１次元 SW(delta)

図 4.１次元 TDSE 比較

図 5.差分法による 2D-TDSE

図 3.２次元 HW(gaussian)

図 6.SW による 2D-TDSE

1P-117

データマイニングによる塩水濃度と水溶液の動的構造に関する研究(北陸先端大院・知識) ○杉山 歩, 中村 大城, 山口 和宏, 水上 卓,
Dam Hieu Chi, Ho Tu Bao	 	 【 序 】塩水は自然界に大量に存在する私たちのもっとも身近な水溶液の一つである。例えば、海水は地球上における水分の97%を占めており、また医療や生体化学の面でも重要な溶液となっている。また、工業的な側面からは塩水は摺動部に腐食摩耗が発生し、機器の性能や寿命に悪影響を及ぼすことから腐食摩耗と電気化学的腐食および両者の相互作用対する研究が進められている。	
 	  現在、塩水の研究手法は様々方法が用いられているが、その一つとして分子シミュレーシ
ョンが挙げられる。分子シミュレーションは方程式を解くことでコンピュータ上に物理特性その特性を示し、解明する手法である。近年の計算機の急速な発達により大規模かつ長時間のシミュレーションを実行することかが可能となり、水様液系に至っては数万〜の水分子を配した系をマイクロ秒オーダーでの計算も可能となっている。	
 	
  一方、近年これらの巨大データから統計的手法により帰納的に法則を見つけ出す、データ駆動型アプローチによる研究が複雑系分野を中心に大きなトレンドとなりつつある。そこで、本研究では巨大データの有効的に使用するため帰納法的アプローチにより新しい知識を獲得する手法の開発し、塩水に対するシミュレーションデータから分子ダイナミクスに関する法則を見つけ出すアプローチを採用する。	
 	
  巨大データの帰納法的アプローチとしてデータマイニングが挙げられる。データマイニン
グは主にマーケティング等の統計処理で広く用いられており、特にWeb上の巨大な情報が複雑なネットワークの解析手法として成果を上げている。例えば、検索エンジンではビッグデー
タの解析手法としてデータマイニングを用いて設計されている。これは計算機に多量データから学習させることにより、	
 内在する法則性を導き出す。	
 	
  このデータマイニングを用いることでシミュレーションにより得られる巨大データを解析することで、従来と異なるアプローチで塩水特性の解明に取り組む。	
 	
 【 計 算 モ デ ル 】本計算はAmber10(amber force field03)を使用し、TIP3P, TIP4PEW水分子モデルを配し、塩水濃度をTable 1に示したNa+, Cl-イオンの数によりコントロールし、NaCl
水様液系及びbulk水系に対し分子動力学計算を行い、ターゲットデータを作成した。解析には先行研究で用いたデータマイニングによるタンパク質水様液系の構造解析手法を採用した[1]。	
 データマイニングには水分子の座標データの時間経過による変化を追うことで、水分子の振る舞いを定量化し、解析をおこなった。解析方法は以下のプロセスに依る。	
 	
 
	
 	
 	
 1.	
 水分子の軌道をステップ数で等分する。 	
 Total atom

	
 	
 	
 2.	
 等分した軌道の重心を求める。	
 
	
 	
 	
 3.	
 重心間距離を求める。	
 等分を細分化する事で一つの分子が	
 64	
 等分され、重心間距離が求めるまで繰り返す。	
 ここでそれぞれの分割ごとに得られた重心間距離に対して次の値を求める。	
 
•	
 平均値	
 
•	
 標準偏差	
 	 •	 最大値 	
 •	
 最小値	
 上記の計23個の値を用いて水一個の分子挙動を一つの	
 23	
 次元のベクトルとして表現し、そのベクトルを主成分分析に正規化したものをプロットする。	
 【 結 果 】上記の手法によって計算した各塩分濃度の水様液の全水分子の第一及び第二主成分の解析結果を図１に示す。詳細は当日報告する。
	
  	
 

	
  	
  	
  	
  	
 

図１：Bulk 水と NaCl 溶液中の水分子のダイナミクス分類結果

References [1] Hieu-Chi Dam, Tu Bao Ho, and Ayumu Sugiyama. 	
  IJCAI, pp. 2424–2429, 2011.

1P-118

Cr(II)二核錯体の hybrid DFT 計算におけるHF 項の混合割合と汎関数の評価（阪大院理 1・理研 2）北河 康隆 1・松井 亨 1・宋 鍾元 1・畑ヶ 宇宙 1・安田 奈都美 1・
川上 貴資 1・山中 秀介 1・奥村 光隆 1

HF ratio and functional set for hybrid DFT
calculation of di-Cr(II) complex(Osaka Univ.1, RIKEN) Yasutaka Kitagawa, Toru Matsui, Jong-Won Song, Hiroshi Hatake,Natsumi Yasuda, Takashi Kawakami, Shusuke Yamanaka, Mitsutaka Okumura【序】 近年の手法および計算機の進歩は、実在する物質の様々な物性量の定量的な量子化学計
算を可能とした。主として多核遷移金属錯体などの擬縮退した分子系における磁気的相互作用は、分子の構造・電子状態を鋭敏に反映し、且つ実験でも観測可能な事から、その定量的計算は磁性物質の研究のみならず、反応中間体や化学結合の安定性を見積もる上でも大変重要となる。一般
的に、分子間および分子内の磁気的相互作用は、ハイゼンベルグハミルトニアンの交換積分（J）値で議論される。

Hˆ = −2∑ J ab Sˆ a ⋅ Sˆ b
巨大な多核遷移金属錯体の金属イオン間に働く磁気的相互作用の計算には、post Hartree-Fock(HF) 法などは計算機コストの問題から適用が難しく、密度汎関数（DFT）法を Broken-symmetry 法と€ともに使用する場合が多い。従って、BS-DFT 法による J 値の定量的算出が課題となる。当研究グル
ープでは以前より、HF 交換項を含んだ hybrid DFT 法が J 値の計算には有効で、その混合割合が重要である事は指摘してきたが[1]、様々な汎関数での系統的な考察は行ってこなかった。そこで、本研究では、Cr(II)２核錯体（図１）の分子内磁気的相互作用に着目し、現在提案されている様々な汎関数での HF 交換項の割合と J 値を比較検討した。

【計算】 錯体は、図１に示したクロム(II)２核アセテート錯体：
Cr2(O2CCH3)4(OH2)2 (1)に着目した。擬縮退した電子状態の計
算には BS 法を用い、J 値は、スピン分極型の１重項(BS)と９重項(HS)とのエネルギー差から山口式

図１ Cr2(O2CCH3)4(OH2)2 錯体

J=

E BS − E HSHSBS− Sˆ 2Sˆ 2を用いて算出した。基底関数は Cr(II)には MIDI+pd、その他には 6-31++G**を使用した。

€【結果】 まず X 線構造に基づき、様々な汎関数で J 値の計算を行った。得られた結果を HF 交換項の割合と J 値の関係でまとめたものが図２である。このように、様々な交換相関汎関数が提案されているが、大まかには、J 値は HF 交換項の割合に比例する事が明らかになった。続いて、交換汎関数
として B88 を、相関汎関数として LYP を用い、HF 交換項と B88 交換汎関数の比率を、Adiabatic
connection（AC）法と Long-range Correction (LC) 法において混合係数を変化させながら J 値を求めた。その結果は、図３のようになり、X 線構造を用いた場合に、実験で得られた J 値を再現する HF
交換項の係数は、AC 法で 0.566、LC 法で 0.720 であった。

HF 法
Exp.=–490cm–1

図２ X 線構造で計算された J 値

図３ HF 混合係数と J 値

しかしながら、そのパラメータを用い、Cr 原子の座標をスピン射影構造最適化法で最適化すると、
Cr-Cr 距離が AC 法で 2.493Å、LC 法で 2.501 Å と実験値 2.362 Å よりも過大評価した。つまり、X線構造で J 値に対して最適化した係数は、構造最適化には使用できない。最適化構造と J 値との両方を満足する混合係数を見いだす為に、AC 法で混合係数を変えながら、最適化構造での J 値をプロットしたところ、実験値ともっとも一致する係数は、構造では 0.37、J 値では 0.47 となった。つまり、
B88 と LYP を使用する限りにおいて、J 値と最適化構造との両者を満たす混合係数は存在せず、汎関数自身を考慮する必要がある事が分かった。当日は、他の汎関数での結果も合わせて報告する。

Reference[1] Y. Kitagawa et al. Int. J. Quant. Chem., 2001, 84, 592-600.

1P119多成分分子軌道法のプログラム開発と応用計算(九大稲盛セ)○石元 孝佳・古山 通久【序論】
水素結合系やプロトン(水素)移動反応など、多くの実験結果から原子核の量子力学的取り扱いの重要性が指摘されている。中でも水素(H)と重水素(D)の違いが引き起こす同位体効果は、水素結合長(Ubbelohde 効果)や水素結合系誘電体の構造相転移など多くの分野で知られており、この H/D
同位体効果は化学反応過程の解析などにも積極的に利用されている。特に水素結合が重要な役割を担っているたんぱく質などの生体内分子では、その機能解明に向けて水素・重水素置換に伴う構造の安定性や反応性に関する詳細な実験的な解析が盛んに研究されているが[1]、溶媒としてのH2O と H2O の沸点や融点といったマクロ物性の違いなど明確に示されていないのが現状である。
我々は H/D 同位体効果を記述するために原子核(プロトン)の量子効果を露に考慮した多成分分子軌道(MC_MO: multi-component molecular orbital)法を開発している[2]。すでに我々は、この
MC_MO 法が H/D 同位体効果の解析に有効であることを示してきたが、実在系を効率よく解析するためにはさらなるプログラムの開発が必要である。そこで本研究では、密度汎関数理論(DFT)や摂動理論に基づく MP2 レベルでの効率的な電子相関の取り込みを実現するためのプログラム開発を行った。また、溶媒効果の取り扱いを実現するために、MC_MO 法に PCM 法を実装した。

【方法】
量子化学計算プログラムである GAMESS をもとに、DFT と MP2 に対応した MC_MO プログラムの実装を行った。今回のプログラム開発では電子相関のみを考慮し、種々の DFT 手法の取り扱いを可能としたプログラム開発に主眼を置いた。また、MC_MO 計算における溶媒効果の取り込みを実現するために、本研究では PCM 法に着目した。電子の基底関数には 6-31G*、プロトンの基底関数には[1s]、[1s1p]、[1s1p1d]GTF を用いて、MC_MO 計算を実行した。

【結果】
HF、B3LYP、MP2 レベルの MC_MO 計算で得られたエネルギーを Table 1 に示す。この結果はプロトンの基底関数として[1s]GTF を設定して得られた値である。ここでは参考のため通常の MO 計算で得られた値も示した。HF に比べて B3LYP や MP2 のエネルギーは低くなっている。また MO 計算と比較すると電子相関による安定化エネルギーが同程度であった。HF レベルで PCM を考慮した
MC_MO 計算では約 0.01 Hartree ほどエネルギーが低くなった。PCM を考慮した MO 計算でも同程度の安定化がみられることから MC_MO 計算でも PCM による溶媒効果が適切に考慮されている

と考えられる。PCM を考慮した B3LYP 法は現在開発中である。

Table 1 Total energy of water molecule. Unit is Hartree.
HF
B3LYP
MP2
MC_MO
Gas
-75.918073-76.270212-76.106501
PCM-75.928544
MO
Gas-76.010747-76.371821-76.196848
PCM-76.021033-76.381127
Fig. 1 には水分子の最安定構造を示す。B3LYP を除いた全ての計算で MO より MC_MO 計算で得られた O-H 間距離が長くなっており、非調和性が考慮されていることがわかる。B3LYP レベルの
MC_MO 計算で O-H 距離が短くなる原因については現在解析中である。
0.9601

0.9472

HF

0.95120.95850.9497

HF(PCM)

0.9691

B3LYP

0.97730.9689MP2

Figure 1 Optimized structure of water. Bond distances obtained
MC_MO and MO are shown left and right, respectively.
沸点や融点などの H2O と D2O 間で見られる水物性の違いを解析するために Fig. 2 に示したような水クラスターを取り上げ、水素結合強度や構造の違いを検討中である。結果の詳細については当日報告する。

謝辞
本研究の一部は平成 24 年度「産学官地域連携による水素社会実証研
究」の助成により行われた。関係各位に感謝する。

【参考文献】

Figure 2 Model structure of (H2O)20 and (D2O)20.

[1] Z. Shi, C. A. Olson, N. R. Kallenbach, and T. R. Sosnick, J. Am. Chem. Soc., 124, 13994 (2002).
[2] T. Ishimoto, M. Tachikawa, and U. Nagashima, Int. J. Quantum Chem., 109, 2677 (2009).

1P120プロトン移動反応の量子的経路について（原子力機構）○志賀 基之【序】一昨年の分子科学討論会では，量子統計力学に基づいた最小自由エネルギー経路（ここでは“量子的 IRC”と呼ぶことにする）を求めるための新しい第一原理シミュレーション手法を提案した[1]．これは，古典系の反応経路を特徴づけるものとして広く受け入れられている固有反応座標（intrinsic reaction coordinate, IRC [2]）のアイデアを，そのまま量子系に一般化したものである[3]．ただし，量子的 IRC の計算は，通常の古典的 IRC のように簡単に行く
わけではなく，第一原理経路積分法[4]とストリング法[5]を組み合わせて反応経路を最適化す
るので，計算量が大規模となる．本研究では，簡単なプロトン移動反応系を対象として，計算を実行し，反応経路の量子効果，温度効果，同位体効果を系統的に考察する．
【理論】N 個の原子の重心（セントロイドという[6]）の質量加重座標 q1 ,, q N  を変数とし，その座標空間における自由エネルギー面 Aq1 ,, q N  を次のように設定する [7,8]．

Aq1 ,, q N    1 log  q1 ,, q N 
ここに，  q1 ,, q N  は座標 q1 ,, q N  に関する確率分布である．Trotter 数を P とした経路
積分表示において，有効ポテンシャルは [9]
P
N m
1

2
Veff r1,1 rN ,P     i  P2 ri ,s  ri ,s 1   V r1,s ,rN ,s Ps 1  i 1 2
である．ここに， ri, s は虚時間スライス s における原子 i の核座標である．これを用いて
N

1 P

mi ri ,s  q i  P s 1
 PI

 q1 ,, q N   lim   
P 

i 1

3P
 N

N
P
2


m
P

1
i

 exp  V r     1
 Z lim  
d
r
d
r

d
r
mi ri ,s  q i 
i ,1 
i,2
i ,P
eff

2



P  
i 1 2π   i 1  P s 1 と表される． mi は原子の質量， Z は分配関数，  P  P  ，   1 kT である．
自由エネルギー面 Aq1 ,, q N  において鞍点から極小点へ至る最急勾配の経路を量子的IRC と定義する．一般に Aq1 ,, q N  面の形には量子効果と温度効果がともに含まれているため，温度条件や同位体置換により経路が変わる．ただし，古典的極限   0 または高温極限   0 では， Aq1 ,, q N  は余分な定数を除いてポテンシャル面 V q1 ,, q N  に等しくなるため，量子的 IRC は通常の古典的 IRC に一致する．
量子的 IRC は探索するのに，次のような最適化アルゴリズムを採用した．すなわち，経路



 



に沿ったセントロイド変数の M 個の組（イメージという[5]） q11 ,, q1N ,, q1M ,, q M
N を発生させ，①それぞれにかかる力の平均値を第一原理経路積分法で計算し，②力をもとにストリング法によってイメージの位置を更新する，というものである．イメージの位置が動かなくなるまで，①と②を交互に繰り返す．
【結果】まず，分子内プロトン移動の例として，NH3 の反転反応を対象とした．図１は，NH3の反転の実際の計算のスナップショットである．イメージ数を M=41 とし，古典的 IRC 上を初期配置とした．それぞれ Trotter 数 P=16 とした ab initio 経路積分ハイブリッドモンテカル
ロシミュレーション[10]を 500 ステップ行うごとに，ストリング法によりイメージの位置を更新した．イメージの更新は 40 回行ったが，後半の 20 回はほとんど動いていない．図２は最適化された反応経路における自由エネルギープロファイルである．温度上昇や重水素置換にともない量子効果が小さくなり，古典のプロファイルに近づく様子がわかる．

図 1：スナップショット
6

A [kcal/mol]

A [kcal/mol]

64NH3 (P=1)
ND3 (P=1)
NH3 (P=16)
ND3 (P=16)
VT S-VEQ

2

0
0

10

20
30
polymer number

40

4NH3 (classical)
NH3 (300 K)
NH3 (900 K)

2

0
0

10

20
30
polymer number40

図 2：自由エネルギープロファイル
次に，分子間プロトン移動の例として，N2H5－の反応のシミュレーションを行った．計算条
件は上述のものとほぼ同様である．この場合は，量子効果によって最適な経路が短くなり，
自由エネルギー反応障壁が劇的に下がる．詳しい結果は当日報告する予定である．
[1] 志賀基之, 藤崎弘士, 分子科学討論会 2010, 4E06.
(1981).

[2] K. Fukui, Acc. Chem. Res. 14, 363[3] M. Shiga, H. Fujisaki, J. Chem. Phys. 136, 184103 (2012).

[4] M. Shiga, M. Tachikawa,

S Miura, J. Chem. Phys. 115, 9149 (2001). [5] D. Marx, M. Parrinello, J. Chem. Phys. 104, 4077
(1996).

[5] W. E, W. Ren, E. Vanden-Ejinden, Phys. Rev. B 66, 052301 (2002); J. Chem. Phys. 126,

106103 (2007).

[6] R. P. Feynman, A. R. Hibbs, Quantum Mechanics and Path Integrals(McGraw-Hill, New York, 1965).

[7] M. J. Gillan, Phys. Rev. Lett. 58, 563 (1987).

Jonsson, Phys. Rev. Lett. 72, 1124 (1994).
(1981).

[8] G. Mills, H.

[9] D. Chandler, P. G. Wolynes, J. Chem. Phys. 74, 4078

[10] K. Suzuki, M. Tachikawa, M. Shiga, J. Chem. Phys. 132, 144108 (2010)

1P121

自由エネルギー面上の反応経路の探索
RIMED 】 麻田俊雄
【 阪府大院・理、RIMED

Approach to explor
xploree reaction pathways on
the free energy surface(Osaka Prefecture Univ.) Toshio Asada[序] 溶液中や生体内で生じる多くの化学反応において、溶媒効果が重要な役割を担っている。
これまでに、分子集合体における化学反応を取り扱うことが可能な QM/MM 法、自由エネル
ギー勾配(FEG1)法および最小エネルギー経路を見つける方法のひとつである nudged elasticband(NEB)法を用いることで、自由エネルギー面上における反応経路を最適化する方法を提案してきた 2。一方、取り扱う分子が複雑になるにつれて単純最適化で得られる反応経

TS F//

路は初期経路に依存したローカルな局所安F-

定経路となる可能性が高くなる(図 1)。そこで初期経路に依存しない反応経路を得る目
的で、自由エネルギー面上の MD シミュレーションを活用することを計画した。自由エ

reactant

ネルギー面における拡張 MD シミュレーシ

product

ョンを行い、
自由エネルギー面上における広図 1 自由エネルギー面における反応経路最範囲な反応経路探索法を検討したので報告

適化と最適化経路のローカルトラップ問題

する。
[計算方法]
自由エネルギー面上の広範囲な構造探索が目的であるため、自由エネルギーを通常のエネ
ルギーと置き換えることでハミルトニアンを構築した。ハミルトニアン自体には特に物理的意味は持たせていない。
p i2+ G (r1 ,⋯ , rN )i =1 2miN

H ' (p, r ) = ∑(1)

一方、QM/MM 法において QM 領域の FEG は∆GiQM (r1 ,⋯ , rN ) =∂E (r1 ,⋯ , rN )
∂riQM
MM

(2)

で表され、MM 領域のアンサンブル生成が必要となる。そこで、 FEG 計算の高速化のため
MM 領域の時間発展に起因する QM 領域の分極は、10fsec ごとに chelpG 電荷を更新する方法(interval sampling 法)により近似した。QM 計算には RHF/6-31G(d)レベルを、また MM

領域には Amber99 力場を用いた。計算対象とする分子は、すでに多くの研究がなされ水中の自由エネルギーマップが報告されているアラニンジペプチドとし、自由度とし てRamachandran plot に用いられる 2 つの二面角に着目した(図 2)。
φ = 180 .0
ψ = 180 .0a)b)φ = − 131 . 8
ψ = − 11 .1

C
C'

N
φ

CA

N'
Ψ

図 2 アラニンジペプチド分子の a) 二面角と b) 水中の自由エネルギー最適化構造。
ここで、 φ = ∠C − C A − N − C ' および φ = ∠N − C A − C − N 'である。

図 3 に full QM/MM MD シミュレーションと、interval sampling を実行した[結果と考察]場合の FEG の x,y,z 座標成分の収束性を示した。定量的には再現しないものの重要な成分についての FEG は傾向が両者で一致していることから、interval sampling の FEG を用いる
15

15

10

10
force / kcal.mol -1A

force / kcal .mol -1A

ことで最適化において計算時間を大幅に節約することが可能である。

5
0
-5
-10

0
-5
-10

-15
0

1000

2000

3000

time / fsec

図3

5

4000

5000

-15
0

1000

2000
3000
time / fsec

4000

5000

5kcal/molA 以上の FEG 座標成分についての収束性。a) は full QM/MM

MD シミュレーションによる結果。b) interval sampling(10fsec)による結果。
10psec の平均力を算出するために、full QM/MM 法で 4 時間程度要するところ、interval
sampling 法では 4 時間程度で計算が完了する。最適化構造を両者で比較したところ、いずれの二面角においても 1 度以内の誤差で一致した。このことから、図 3 にあげた FEG の誤差は収束する構造にほとんど影響しない程度であると結論できる。
さらに、自由エネルギー面上における拡張 MD シミュレーションを実行し、自由エネルギ
ー曲面上のダイナミクスについて解析を行った結果については当日に発表する。
[参考文献] 1. M.Nagaoka, N.Okuyama-Yoshida, and T.Yamabe, J.Phys.Chem.A. 102, 8202 (1998).
2. N.Takenaka, Y.Koyano, Y.Kitamura, T.Asada, and M.Nagaoka, Theo.Chem.Acc., 130, 215-226(2011). 櫻井耕司、麻田俊雄、小関史朗、分子科学討論会 2012 ポスター 4P113.

1P122補体たんぱく質 C3 が抗原を捕捉する際の結合機構に関する理論研究(神奈川大・理) 笹本千怜, 志村亮, 松原世明

Theoretical study on the process of bonding of the complement protein C3
with the antigen(Kanagawa Univ.) Chisato Sasamoto, Ryo Shimura, Toshiaki Matsubara【緒言】 補体は生体内の免疫システムで重要な役割を果たすタンパク質であり、20 種ほど存在する。免疫システムには、自然免疫と獲得免疫があり、抗原が体内に侵入すると、前者に属す
るマクロファージや後者に属するキラーT 細胞などが抗原を攻撃し、外敵から体を守っている。
特に、マクロファージは、抗原を捕食し消化してしまう。その際、マクロファージは、抗原と結合した補体を認識することで食作用を促進する。食作用を促進させる補体のこの作用をオプソニン化という。オプソニン化は、補体の重要な役割の一つである。本研究は、このオプソニン化に関与する補体 C3 に着目し、C3 が何故抗原をうまく捕捉できるのか明らかにすることを目的とし、量子化学計算により解析を行った。
【計算方法】 補体は全て、活性化の過程を経て機能を発現する。オプソニン化に関与する補体
C3 も、活性化され C3b となり機能を発現する 1)。その際、活性サイトを含む TED と呼ばれるドメインが移動するとともに結合部位が出現する。TED ドメインの活性サイトでは、Cys988 とGln991 が接近し反応サイトとなるチオエステル-S-CO-を形成する（図 1）。従来の反応機構では、そこに、遠くにいた His1104 がやって来て、His1104 のイミダゾール基と Gln991 のアシル基との間で相互作用する。その後、S-C 結合が開裂し、抗原の OH 基を認識して抗原を捕捉すると考えられている 2)。そこで、本研究では、この反応機構に含まれる中間体と遷移状態の構造最適化を行ない、反応のエネルギープロフィールを求めた。モデル分子として、活性サイトについては
Cys-Gln-His、抗原については、実験でもモデル分子として採用されているアルコール（メタノール）を用いた。構造最適化は、密度汎関数法(B3LYP)により、6-311++G(d,p)基底系を用いて行った。平衡構造および遷移状態は、振動数を計算し確認した。反応座標は、IRC (a) 従来の反応機構
His1104
計算により確認した。さらに aug-cc- Cys988
+
+
−
S C N
SH C NS C N
pVTZ 基底系を用い、エネルギーを計
SH C N
O
RO O
−
O
O
R
O
算した。熱力学的パラメータの計算は、
ROH
Gln911(抗原)
振動数を用い、298.15 K で行った。
全ての計算は、GAUSSIAN03 プログ (b) 新規反応機構
N
N
N
N
N
N
ラムを用いて行った。
H1 H2
H H
H H
TS2n
TS1n
【結果と考察】 従来の反応機構に含
O1
OO
S Cys
Gln
S Cys
S Cys
Gln
Gln
まれる遷移状態および平衡構造の構造
O
O
O H
H
R
H
R
R
最適化を行った結果、この反応機構は、
N
N
N
(抗原)
4 中心の遷移状態（図 2）を経由する
His
His
His
N1
N
N
Pn
In
Rn
H
H
H3
１段階の反応であることが分かった。
O - O2
O-O
O -O
Asp
Asp
Asp
アミノ酸の His は、Gln のカルボニル炭素と弱く相互作用しているものの、
図 1． 補体 C3b が活性サイトで抗原を捕捉する際の反応にはあまり関与しないことが分か
従来の反応機構と新規反応機構
った。したがって、C-S 結合は His の2.545 O
S
C
1.688
1.801
H
O
Unit in Å 1.095

0.0
Rn

5.0

36.6
TS1o

6.8

4.0
TS1n
新規反応機構4.0

従来の反応機構

-2.4
In

TS2n

-6.5 Pn-12.7

Po

Reaction coordinate

図 2．補体 C3b が活性サイトで抗原を捕捉する際の従来の反応機構と新規反応機構のエネルギープロフィー
ルおよび前者の遷移状態構造

Distance (Å)

Energy (kcal/mol)

影響をほとんど受けることなく、C-S 炭素はメタノールの O-H 酸素の求核攻撃を受ける。遷移状
態では、O-H 酸素と C-S 炭素の距離は既に短くなっており、結合がほぼ生成していることが分かる。遷移状態の後、O-H 水素はプロトンとして硫黄に引き取られ、反応は完了する。この反応は、わずかな発熱反応であり、活性化エネルギーは、36.6 kcal/mol と非常に大きいことが分かった（図2）。また、His の影響によってエネルギー障壁はわずか 1 kcal/mol しか低下していないことが分かった。
このような大きなエネルギー障壁を越えて反応が進行しているとは考えにくい。そこで我々は、図 1(b)に示す従来の反応機構とは異なる新たな反応機構を考えた。これは、セリンプロテアーゼによるたんぱく質分解反応の反応機構と同様の反応機構である。この場合、His は反応を仲介する役割を果たす。His がメタノールの O-H 水素をプロトンとして引き取ると同時に O-H 酸
素が C=O 炭素に付加し、中間体 In を生成する。この際、O-H 酸素の負の電荷は C=O 酸素に移動する。その後、C-S 結合が開裂するとともに、His は引き取ったプロトンを硫黄に引き渡す。
この反応機構の反応物 Rn、中間体 In、生成物 Pn および最初の過程 Rn→In の遷移状態 TS1n、その後の過程 In→Pn の遷移状態 TS2n を構造最適化しエネルギープロフィールを求めた。図 2に示すように、２つのエネルギー障壁はともに小さくエネルギープロフィールは滑らかである。
その理由の一つは、中間体 In がエネルギー的に大きく安定化されるためである。In の C=O 酸素の負の電荷は、周囲のアミノ酸（今回は水分子をモデル分子として使用）の正の部分電荷との相互作用によって安定化される。実際、図 3 に示すように、In は正の部分電荷を持つ２つの H（今回は水分子）と水素結合を形成する。また、基質と His との間のプロトンの受け渡しでは、Aspが重要な役割を果たしている（図 3 では Asp の周囲のアミノ酸の効果は考慮していない）。これ
らのアミノ酸の効果は、これまでセリンプロテアーゼの反応機構で示唆されてきたことであり、我々は他の類似の酵素反応においても量子化学計算によって確認してきた 3) 。また、メタノー
ルの O-H 酸素の求核性は、His との相互作用によって高くなっていることも分かった。このように、補体 C3b の抗原との結合反応は、従来の機構ではなく、我々が提案する新たな機構で進行すると考えられる。当日は、エステル結合がチオエステル結合であることの重要性、つまり S 原子をもつ Cys の寄与の重要性、さらにその他の周囲のアミノ酸の効果なども議論する予定である。

H1-O1

3.0

H2-O1

2.0
1.0
0.0

N1-H3

H3-O2

Rn TS1n In TS2n Pn

図 3．新規反応機構の反応座標に対する構造パラメータの変化【参考文献】
1) A. A. Ajees, K. Gunasekaran, J. E. Volanakis, S. V. L. Narayana, G. J. Kotwal, and H. M. K. Murthy,
Nature, 444, 221-225 (2006).
2) S. K. A. Law and A. W. Dodds, Protein Science, 6, 263-274 (1997).
3) Y. Sakae, T. Matsubara, M. Aida, H. Kondo, K. Masaki, and H. Iefuji, Bull. Chem. Soc. Jpn., 82,
338–346 (2009).

1P123
遷移金属原子を含んだ巨大系を取り扱うための IMiC-MO 法の拡張
(岐阜大工)

宇田川太郎, 酒井章吾

Application of IMiC-MO method to transition metal complexes(Gifu University) Taro Udagawa, Shogo Sakai【序論】

遷移金属原子は生体内において生体高分子と錯体を形成し、重要な生理機能を担っていることが知られている。また近年では、例えばハロゲン架橋錯体のような、多彩な興味深い物
性を示す無機錯体にも注目が集まっており、様々な分野で盛んに研究されている。これら遷移金
属錯体の興味深い物性の多くは、d/f 電子の存在により引き起こされるが、それらは同時に理論化
学計算による取り扱いを複雑にしている。このような系では、計算結果は初期分子軌道に強く依存する。また、生体内高分子錯体や集積型複核錯体などのような巨大系に対しては、従来の分子軌道法による直接的な取り扱いは現実的ではない。
近年我々は、系を分割することで系のエネルギーおよび各原子にかかる力を効率良く求めることが可能な Integrated Multicenter MO (IMiC-MO)法を開発してきた[1]。この手法では、全系を目的
分子、隣接分子、環境分子の 3 つの領域に分割することで、効率良く巨大系を取り扱う。近年では様々な分割的手法[2]が提案されているが、それらの手法を遷移金属錯体に適用した例は少なく、特に電子状態が複雑となる少数の d 電子を有する系へと適用した報告は無い。本研究ではまず岩田により開発された球状電荷解析法[3]により、古典的酸化数と従来の電荷解析法、および実際の遷移金属錯体周囲の電子分布について明らかとし、金属錯体を分割法により取り扱う正当性を示した。そして種々の単核および複核錯体において d 電子の様々な占有状態を IMiC-MO 法で取り扱い、複雑な電子状態を持つ遷移金属錯体を解析した。

【計算方法】

遷移金属錯体のモデルとして、Ti, V, Cr, Mn, Pt の単核および複核錯体を用いた。

中心金属原子の酸化数は、Ti は II~IV, V は II および III, Cr は II, Mn は I, Pt は II とした。構造はそれぞれの錯体の最も安定な電子状態において最適化した。全ての計算は(U)B3LYP 法で行い、遷移金属錯体に LanL2DZ、その他の原子には 6-31G を用いた。従来の(U)B3LYP 法により最適化した構造に対して IMiC-MO 法では、同じ計算方法/基底関数を用いてエネルギーを評価した。

【結果・考察】図 1 に球状電荷解析法により求めた、種々の錯体における電子の平均密度差Δρ0(R)を示した。ここで、Δρ0(R)は遷移金属原子を中心とした半径 R の球表面の電子密度の、錯体形成による変化を表している。図より、遷移金属原子が異なっても酸化数の同じ錯体では良く似た Δρ0(R)を示すことがわかった。これは中心金属原子が正に帯電することで、配位子の電子をより強く引きつけるためであると考えられる。また、IMiC-MO 法の手順では、遷移金属錯体の金属
原子単体を目的分子として取り扱うため、配位子との電子移動量に相当する Δρ0(R)が分割誤差に

直接影響すると考えられる。Ti(IV)錯体では、

0.18

Δρ0(R)の最大値は他の錯体に比べ非常に大き0.16

Ti(III)Cl4
Ti(III)Cl4

0.14

い の に 対 し て 、 Ti(II) 錯 体 で は 全 域 に 渡 りΔρ0(R)は非常に小さい。これは Ti(II)原子周囲の電子状態が錯体形成により殆ど変化しておTi(II)Cl4
Ti(II)Cl4Δρ0

Ti(IV)Cl4
Ti(IV)Cl4

0.12

Ti(II)(H2O)6
Ti(II)(H2O)6

0.10

Ti(IV)(H2O)6
Ti(IV)(H2O)6

Ti(III)(H2O)6
Ti(III)(H2O)6
V(II)(H2O)6
V(II)(H2O)6

0.08

らず、+2 の Ti(II)原子様状態で存在しているこ

0.06とを意味している。これら錯体中の Ti(II)原子0.04

V(III)(H2O)6V(III)(H2O)6Cr(III)(H2O)6Cr(III)(H2O)60.02

上の電荷は Mulliken 解析法では約 0.86 であり、

0.00

よく知られているように従来の電荷解析法で-0.02は遷移金属上の電荷を適切に表現できていな

0.0

0.2

0.4

0.6

0.8

1.0

r [Å]

いことがわかる。それよりも、古典的な酸化図1. 種々の遷移金属錯対における電子の平均密度差数から考えたモデルの方が非常に単純ながら
静電場に対する近似として適していると考えられる。
次に、Ti, V, Cr の II および III 価の六水和物に対して、種々の電子配置を従来の(U)B3LYP 法お
よび IMiC-(U)B3LYP 法により計算した。
IMiC-MO 法において目的の電子状態を取り扱うには、初期軌道の与え方が重要である。本研究では、分割した最小単位の系の密度行列から作成した各分%&'()*+,-./0123456789:;<9=>4?7@ABCDEFGH.
/I+JK@ABCD(LMNOPQRSPTABCDEFBG@HB?I0!"#$!@!,JHF-)*+,# D!"#$%&'( )*+,-!"#$%&'(!/ 01# 2 !/ 13 # 2
456%
7577
7577
.+!))#%.2!/ 30 # 2!/ 01# 2
4568
2592
258:
!/ 30 # 2!/ 13 # 2
456%
2599
2599
%

割系に対する近似密度行列を基に各初.;!/ 30 # 2!/ 1; # 2

<58<

;958;

;65:<

.2

!/ 13# 2!/ 01# 2!/ 30 # 2

-285;2-285;8-285;;

7577
259%
2598

7577
;572
;577

.+!)))#;

期軌道を作成した。表より、中心金属原子の種類に関係なく、同じ価数の錯=!))#体においては IMiC-(U)B3LYP 法により4>;!/ 30 # 2!/ 01# 2 !/ 13# 2

-:5<8

7577

7577

4

.;

!/ 01# 2 !/ 13 # 2 !/ 3;-0; # 2-2%52;

%%59;

425;6

%

.2

!/ 01# 2 !/ 13 # 2!/ 30 # 2!/ 01# 2!/ 30 # 2!/ 13 # 2-2<599-285;7-28522

7577
;54;
;5:9

7577
;56%
;5<2

@

!/ 30 # 2!/ 01# 2 !/ 13# 2!/ 1; # 2

-;544

7577

7577

.2

!/ 30 # 2!/ 01# 2 !/ 13# ;!/ 30 # ;!/ 01# 2 !/ 13# 2!/ 30 # 2!/ 01# ; !/ 13# 2

7562
754%
759:

;657%
;95<9
;95<9

;;59<
;659%
;6542

>;

!/ 30 # 2!/ 01# 2 !/ 13# 2-265;;--=!)))#

取り扱ったことによるエネルギー誤差
ΔE は Δρ0(R)から期待されたように同程,?!))#度であり、全体的に精度よく全エネル:%

ギーを再現した。各電子配置間の相対エネルギーについても、(U)B3LYP 法に,?!)))#

4

Dよる結果を非常に良い精度で再現した。 !@!,JHF-)*+,#GKG@!!"#$%&'(#G-G@!)*+,-!"#$%&'(#
その他本研究で取り扱った全ての遷移金属錯体で IMiC-(U)B3LYP 法は従来法による結果を精度良く再現し、IMiC-MO 法により複雑な電子状態を持つ遷移金属錯体を取り扱えることが示された。
【References】 [1] S. Morita, S. Sakai, J. Comp. Chem., 22, 1107 (2001), S. Sakai, S. Morita, J. Phys.
Chem. A, 109, 8424 (2005), T. Udagawa, S. Sakai, J. Comp. Chem., 32, 2902 (2011). [2] K. Kitaura, E.
Ikeo, T. Asada, T. Nakao, M. Umebayashi, Chem. Phys. Lett., 313, 701 (1999), M. Kobayashi, H. Nakai,Phys. Chem. Chem. Phys., 14, 7629 (2012). [3] S. Iwata, Chem. Phys. Lett., 69, 305 (1980).

1P124チミン分子における溶液内超高速無輻射失活の QM/MM‒MD シミュレーション（北大院総合化学1，北大院理2） 新井 岳1，中山 哲2，山崎 祥平2，武次 徹也2Excited-state QM/MM‒MD simulation of thymine in aqueous solution(Hokkaido Univ.)

Gaku Arai, Akira Nakayama, Shohei Yamazaki, and Tetsuya Taketsugu【研究背景】DNA 塩基は紫外領域に吸収帯を持つが、高効率な無輻射失活経路を有しており、この性質は DNA 塩基の光安定性として知られている。近年の精力的な研究により、この失活過程はサ
ブピコ秒のオーダーで起こり[1,2]、ポテンシャル曲面の円錐交差(conical intersection: CI)点が重要な役割を果たしていることが分かっている。気相中の単分子塩基の失活過程に関しては実験的にも理論的にも多くの知見が得られてきたが、溶液内や DNA 骨格内での失活過程に関しては、周囲の環境が励起寿命に及ぼす影響等、未だ十分に理解はされていない。本研究では、水溶液中でのチミン分子を対象とし（図 1）、励起状態 QM/MM‒MD シミュレーションを行い、気相中での計算と比較することで、周囲の溶媒分子が失活過程に与える影響を議論す
る。
(a)

図 1

(b)

(a) チミン分子 (b) QM/MM-MD 計算のスナップショット【計算手法】
チミン分子に対して、気相中(QM)と水溶液中(QM/MM)での MD シミュレーションを行った。QM 領域として取り扱うチミン分子の量子化学計算には、基底状態に対しては MP2 法を用い、励起状態に対しては、気相中の計算では CASSCF(2,2)を、水溶液中では CASSCF(4,4)
を参照関数とした CASPT2 法を用いた（それぞれ CASPT2(2,2), CASPT2(4,4)とする）。基底
関数は全て Sapporo‒DZP を用いた。また、MM 領域である水分子には SPC/F の力場を使用した。水溶液中での計算では、一辺 16 Å の立方体セル内にチミン 1 分子と水 122 分子を配置し
た。

【結果】
チミンの光励起緩和過程には、2つの励起状態(ππ*, nπ*)の関与がこれまでに示唆されているが[3]、本発表では明状態であるππ*状態のみを経由した過程に焦点を当てる。これまでに、
C5=C6のねじれを伴う円錐交差点(ππ*(C5=C6)/gs)CIを経由した失活過程が報告されている。
基底状態で温度一定(T = 300 K)のMDシミュレーションを行い、系を平衡化させた後、一定時間間隔で座標と運動量を取り出し、ππ*状態でのQM(QM/MM)‒MDシミュレーションを開始した。MDシミュレーションでは、計算コストの観点から活性空間の小さいCASPT2(2,2)またはCASPT2(4,4)を用いているが、トラジェクトリに沿って一定の間隔でCASPT2(12,9)のエネル
ギー計算を行い、ππ*状態に関してはこれらの小さな活性空間の計算で十分な精度が得られていることを確認した。
代表的なトラジェクトリの結果を図2に示す。気相中では、励起直後はππ*とnπ*状態のエネ
ルギーが近接しており、ππ*状態の運動では、C5=C6結合のねじれを伴い、(ππ*(C5=C6)/gs)CIに到達することが示された（図2(a)）。また、水溶液中の同様のトラジェクトリでは、気相中に比べてππ*状態が安定化し、nπ*状態が不安定化することが分かった。この結果は、溶液中ではππ*のみが失活過程に重要であり、nπ*状態は関与しないことを示唆している（図2(b)）。
さらに、水溶液中の計算では、気相中で見られなかった、(ππ*(C5=C6)/gs)CIとは異なる円錐
交差点に到達するトラジェクトリが得られた（図2(c)）。この円錐交差点はC4=O8が六員環平面から大きく歪んだ構造をしており、水溶液中ではそのエネルギーが気相中に比べて大きく安定化されていることがわかった。以上のように、本研究で水溶液中でのみ起こる新しい失活経路の存在を示すことができた。トラジェクトリ数を増やして統計的な解析を行った結果については当日報告する。

(a)

(b)

―ππ* CASPT2(2,2)―gs
CASPT2(2,2)●nπ* CASPT2(12,9)●ππ*
CASPT2(12,9)●gs
CASPT2(12,9)(c)―ππ* CASPT2(4,4)―gs
CASPT2(4,4)●nπ* CASPT2(12,9)●ππ*
CASPT2(12,9)●gs
CASPT2(12,9)

図 2: 代表的なトラジェクトリーにおけるポテンシャルエネルギーの時間変化。
(a) 気相中（(ππ*(C5=C6)/gs)CI に到達）、(b) 水溶液中（(ππ*(C5=C6)/gs)CI に到達）、
(c) 水溶液中（(ππ*(C4=O8)/gs)CI に到達）。

【参考文献】[1] C. Canuel et al., J. Chem. Phys. 122, 074316 (2005).
[2] T. Gustavsson et al., Chem. Phys. 350, 186–192 (2008).
[3] S. Yamazaki and T. Taketsugu, J. Phys. Chem. A 116, 491–503 (2012).

1P125直線炭素鎖アルコール HC4OH のマイクロ波分光(上智大学 1、静岡大学 2、東京理科大学 3) 國松亜利沙 1、田辺沙織 1、小川智史 1、久世信彦 1、
中根綾 2、岡林利明 2、荒木光典 3

Microwave spectra of the linear carbon-chain alcohol HC4OH(Sophia Univ.1, Shizuoka Univ.2, Tokyo Univ. of Science3) Arisa Kunimatsu,1 Saori
Tanabe,1 Satoshi Ogawa,1 Nobuhiko Kuze,1 Aya Nakane,2 Toshiaki Okabayashi,2Mitsunori Araki3【序】星間空間には二重結合または三重結合をもつ多くの炭素鎖分子が確認されている。また、ヒドロキシル基をもつ分子も多く確認されている。しかし、これら両方の特徴をもつ分子はこれまで報告されていない。我々はこれら両方の特徴をもつ分子、すなわち、ヒドロキ
シル基をもつ直線炭素鎖分子 1,3-butadiyn-1-ol (HC4OH、図 1)のマイクロ波スペクトルを実験室系で観測してきた[1-3]。これまでは熱分解反応を利用して HC4OH を生成してきたが，本研究では新たに diacetylene を用いた方法を試みた。さらにこの分子の天文観測に必要となる、精密な静止周波数の値を得るため、フーリエ変換型マイクロ波分光器による観図 1. 1,3-butadiyn-1-ol (HC4OH)測も行ったので報告する。

【実験】マイクロ波スペクトルの測定は 100 kHz 矩形波シュタルク変調型マイクロ波分光器と、ファブリペロー型共振器を用いたフーリエ変換型マイクロ波分光器を用いて行った。

シュタルク変調型マイクロ波分光器を用いた実験では、脱ガスを行った 2-butyn-1-ol(H3CCCCH2OH)をサンプル管に入れ、電気炉で 900 ℃まで加熱した石英ガラス管に通すことで熱分解反応を行い、反応生成物をそのまま 3 m の X-band の導波管セルに導入した。重水素置換種は、2-butyn-1-ol と重水を気相で混合した後、熱分解する方法で合成した。マイクロ波スペクトルの測定は、周波数領域 12.0～40.0 GHz、試料圧 10～17 Pa、シュタルク電
圧 30～800 V、積算回数 5～16 回の条件下においてフロー状態で行った。

21200

21300

J=5-4
P=10 Pav=1 E=60 V
積算回数：16回

J' Ka Kc －J" Ka Kc =514－413
積算回数：200回

v=2

intensity

intensity

v=0

21400

21500 21284.8

frequency/MHz

図 2. HC4OH のスペクトル（シュタルク変調）

21285.3

frequency/MHz

21285.8

図 3. HC4OH のスペクトル（フーリエ変換）

フーリエ変換型分光器を用いた実験では、化学的に合成した diacetylene を O2 ガスとともに Ar ガスで希釈しパルス放電による 1,3-butadiyn-1-ol の生成を試みた。測定は周波数領域
8~26 GHz、試料背圧 1.0～1.3 気圧(diacetylene:0.75%、O2:0.75%、Ar:98.5%)、積算回数
50～500 回の条件下で行った。

【結果と考察】フーリエ変換分光を用いたノーマル種の測定では、シュタルク変調分光での観測周波数の近傍においてラインを観測できた。特筆すべきことに、これらのラインは O2がなくても、放電を行わなくても観測された。このことは diacetylene の化学的合成の際に、副生成物として 1,3-butadiyn-1-ol が生成していることを意味する。さらに、それぞれの J において、低周波数側にさらに強度の大きいラインが見つかったため、これらを Ka=0 と帰属した。したがって、シュタルク変調分光でこれまで Ka=0 と帰属されていた吸収線は、Ka=1であることがわかった。最終的に、シュタルク変調分光では 50 本、フーリエ変換分光では15 本 の ピ ー ク を 帰 属 で き た 。 こ れ ら の 遷 移 周 波 数 に 対 し て Watson の A-reduced

Hamiltonian を用いて最小二乗解析し、振動基底状態および振動励起状態の回転定数を求めた(表 1)。振動基底状態ではフーリエ変換分光の結果を、振動励起状態ではシュタルク変調分
光の結果を示す。
一方、重水素置換種のスペクトルは、シュタルク変調分光での測定によって、3 つのグループ（それぞれ I,II,III group とする）が存在することが分かった。各グループのピーク間隔と理論計算による B+C の値の比較から、I group と II group はそれぞれ DC4OH、HC4OD

由来のスペクトルだと判断される (表 2,3)。III group のピーク間隔は、HC4OH の B+C の理論計算値と近いが、観測されたピークの遷移周波数が異なっているため、異性体によるピークの可能性があると考えられる。現在、化学合成された 1,3-butadiyn-1-ol を含む diacetyleneの試料を NaOH/NaOD の 1:1 混合溶液を使って重水素化することで，DC4OH、HC4OD の帰属を確定するための準備を進めている。
表 1. HC4OH の回転定数
a)

c)

v =0

v =1

c)

v =2

b)

B /MHz2126.5060(6)C /MHz2124.0602(6)6874002137.55(3)2134.88(5)687400c)b)

A /MHz6874002133.44(2)2130.63(3)b)∆J /MHz0.003407(8)0.003407

0.003407

a)
a)

b)

∆JK /MHz-2.8093(3)-2.8093-2.8093

a)
a)

c)
a)
a)
6874002141.29(1)2139.68(1)v =30.003407-2.8093a) A 、v =1～3の∆J 、∆JK は固定b) ( )内は1σc) ν =0はフーリエ変換型、ν =1～3はシュタルク変調型の結果を用いた

表 2. 重水素置換体の回転定数(理論値)

表 3. 重水素置換体の回転定数(実験値)

DC4OH

I group

HC4OD

B /MHz1977

C /MHz1972

B +C/ MHz
3949

2027

2015

4042

(計算レベル：MP3/cc-pVTZ)

II group

III group

B +C /MHz
3989

4086

4219

【参考文献】[1] Araki & Kuze, ApJ, 680, L93 (2008)[2] 久世、胡木、荒木、分子科学討論会、名古屋、4P019 (2009)[3] 久世、田邉、小川、胡木、荒木、第 26 回化学反応討論会、広島、1P15 (2010)

1P-126
一重項酸素とβ-カロテンの反応についての DFT 計算	
 (大阪市大・院理学)小泉智也，田中美和，松下叔夫
DFT Study of the Reactions between Singlet-Oxygen and β-carotene(Osaka City University, Graduate School of Science)
Tomoya Koizumi, Miwa Tanaka, Toshio Matsushita【序】カロテノイドは，光合成における重要な色素であることに加え，一重項酸素に代表される活性酸素を除去する能力を有している．この保護機能は，遊離カロテノイドについても見られる．この性質は，カロテノイドのラジカル除去能力に依存しており，一重項酸素を最も有効に失活する化合物のひとつである．失活は chemical quenching と physical quenching の両方が関与すると考えられているが，それらの機構についてはほとんど報告がなされていな
い．特に，chemical quenching が，どのような素反応から構成されているか，共役鎖と酸素の直接付加，あるいはカロテノイドラジカルが関与するのか，実測された酸化物は，一重項酸
素が直接付加したものか，又は一重項酸素からのエネルギー移動により生成した三重項カロ
テンと三重項酸素との反応によるものか，他の競合する反応(ジオキセタンやレチノールへの分解)が失活にどのように影響を与えるのか等，反応の詳細は明らかではない．これらについての知見を得るため分子軌道法計算を行った.	
 【計算方法】計算はすべて B3LYP/6-31G(d)で行い，open singlet については非制限法を用い，
spin contamination は山口らの方法で補正した．β-カロテンのメチル基をすべて水素に置き換えたモデル化合物(以下 H-β-カロテン)とβ-カロテンと酸素の反応を Fig.１のスキームにしたがって計算をおこなった．ここでは，H-β-カロテンについて記載し，β-カロテンの結果
は当日報告する．最も反応性が高いと考えられる共役の末端炭素を 1，カロテンの基本構造における末端の炭素
Fig.1 Chemical quenching reaction scheme for carotene and oxygen
を 3，そしてカロテン
か ら レ チ ノ イ ド へ の 1Carotene + 1O
1Complex
1Biradical
1Dioxetane2
分解に関連している
Further
ISCと思われる共役鎖のreaction
3Complex
3Biradical
3Dioxetane
中心の炭素を 11 とし 1Carotene + 3O2(catalyst)
た(Fig.2) .
Fig.2 Numbering of atoms for H-β-carotene【結果と考察】C1 と O2 の反応におけるポテンシャルエネルギーを Fig.3 に示す．1biradical を生成するためのエネルギー障壁は約 3kJ と非常に小さく， 1biradical は 1complex に比べて約 50kJ 安
定化している．また一重項の 1,2-dioxetane へのエネルギー障壁は約 50kJ となっている．
biradical では一重項と三重項のエネルギー差が約 0.07kJ と非常に小さくなっており，またそれぞれの構造が近いため容易に一重項から三重項への isc による遷移が起こると考えられる．
50kJ 程度のエネルギー障壁をもつ単分子反応の寿命は ms のタイムスケールであり，一方，
一般に biradical の isc における寿命はμs と遥かに小さいことが報告されている．よって
1
biradical からは，一重項の 1,2-dioxetane や 1complex への反応よりも，isc による三重項への遷

移が起こりやすいと考えられる．	
 	
  3biradical から 1β-carotene	
 と O2(3 Σ-g)の 3complex を生成するエネルギー障壁は約 14kJ と低い結果が得られた．この反応は一重項の 1,2-dioxetane や三重項の 1,2-dioxetane を生成する反応と競合するが，エネルギー障壁(それぞれ，50.2kJ と 99.5kJ)を比べると明らかに小さく起こりやすいと考えられる．このことより，biradical で一重項から三重項への遷移が起こり，1β-carotene	
 と O2(3 Σ-g)へと解離するという経路が最も有利であると考えられる.	
 	
  C３ ，C11 との反応の結果を Table	
 1にまとめた．1complex から 1biradical へのエネルギー障
壁は C1 ，C3 ，C11 の順に大きい．biradical における安定性を見ると，C3 は C1 に比べて約 18kJ
不安定化し，C11 についてはさらに 28kJ 不安定となる．3biradical から 3 complex へのエネルギ
ー障壁は C1 ，C3 ，C11 の順に小さくなっており，これは biradical でのエネルギーが高くなるためであると考えられる．これらのことから，もし 10 kJmol-1 程度のエネルギーを反応に使用できるとすると C1~C7 程度まで，つまり共役鎖と O2 との衝突の 2/3 程度が，15 kJmol-1 では
どの共役部分に衝Fig.3 Potential energy profile for the attack on the C1 by 1O2
突 し て も chemical
quenching (1carotene+ 3 O2 への解離)が起こると考えられる．
一 重 項 の 1,2dioxetane，三重項の1,2-dioxetane に お
1complex
け る 安 定 性 もbiradical と同様に，
C1 ，C3 ，C11 の順に不安定になるという結果が得られた．

Table 1 The energy barriers and relative energies in each reaction
Energy barrier
1

1

C1

complex↓1
biradical
3

biradical↓
1
dioxetane
50.2

C3

4.8

C11

15

1)

3

ΔE

2)

3

biradical↓
3complex14.4

biradical↓
3
dioxetane
99.5

50.8

9.6

43.1

1.6

1.3

biradical

1

dioxetane3

dioxetane

44

2.6

60.9

98.2

62

17.5

82.8

120.6

90

33.5

148.5

1) in kJmol-1 2) Relative to 1 β-carotene + O2(3Σ-g), in kJmol-1

1P127電子励起状態における自由エネルギー面の構築: RISM-TDDFT 法(京都大学福井謙一記念研究センター) 嶺澤 範行

Free energy profiles for electronically excited states: RISM-TDDFT approach(Kyoto Univ. FIFC) Noriyuki Minezawa【序】 溶液内分子の自由エネルギーを評価することは、化学反応をはじめとする分子の動的過程を理解する上で重要な知見を与える。特に光化学反応では、分子の電子励起状態を記述することが要求されるため、静的および動的電子相関を考慮した量子化学計算を行う必要がある。現在
QM/MM-MD 法が自由エネルギーを計算する手法として広く用いられている。しかし、精度の良い電子状態計算を長時間のサンプリングと組み合わせて行う場合、計算コストが莫大になる。
RISM-SCF 法は電子状態計算と RISM 積分方程式理論を組み合わせた方法で、溶液内自由エネ
ルギーの計算に幅広く適用されてきた。この方法の利点は、(1)長時間のサンプリングを必要としないため計算コストが小さいこと、(2)変分原理に基づいており自由エネルギーが極小であることである。本研究では、RISM-SCF 法の拡張として、溶液内分子の電子励起状態における自由エネルギ
ー計算を行う手法を考える。そこで、LR-RISM-TDDFT 法を導入する。すなわち、線形応答(Linear
Response)の枠内で、RISM と時間依存密度汎関数法(TDDFT)を組み合わせた方法である。TDDFT
法は計算コストの面で優れ、電子励起状態をある程度定量的に記述できることが知られている。したがって、RISM との組み合わせは、溶液内における電子励起状態の自由エネルギーを効率よく計算するための自然な拡張である。また、LR-RISM-TDDFT 法の解析的な自由エネルギー微分法を導入することで、自由エネルギー面の構築も可能である。
【方法】 LR-RISM-TDDFT 法の導出の詳細は当日に報告する。Solvated Fock 演算子の電子密度依存を考慮して、TDDFT の方程式を修正した。下図に LR-RISM-TDDFT 法の計算の流れをまとめた。線形応答を用いるため、RISM-SCF の計算は、基底状態に対してだけ行うことになる。

【結果】 テストとして、アセトンの n*励起状態を計算した。電子状態計算は B3LYP 汎関数および
cc-pVDZ 基底関数を用いた。溶媒は水(SPC モデル)である。最適化構造を表にまとめた。比較として、気相および誘電体モデルの計算結果も示した。面外角が誘電体モデルよりも若干

図: LR-RISM-TDDFT 法

大きな値を示しており、水素結合のような局所的な溶質‐溶媒間の相互作用の重要性を示唆している。
【考察】 今回導入した LR-RISM-TDDFT 法により、大規模分子の自由エネルギー面の探索が可能となった。しかし問題点として、線形応答の妥当性があげられる。それは右図から明らかなように、溶媒和構造が励起状態の電子分布に基づくものではないことである。これを解決するためには、励起状態の電子密度をあらわに求め、それに対して RISM-SCF 法を適用しなければならない。余裕があれば、この特定の励起状態(state specific)の電子密度を考慮した SS-RISM-TDDFT 法についても言及したい。

表: アセトンの n*励起状態の最適化構造 (TD-B3LYP/cc-pVDZ)
Gas

PCM

RISM

r(CO)/Å

1.316

1.316

1.309

r(CC) /Å

1.508

1.513

1.521

CCC

120.2

119.8

117.8

CCO

114.9

114.6

114.5

O-CCC

32.3

33.7

36.5

Dipole/D

1.85

2.35

2.32

1P128RT 並列法によるメニーコアプロセッサ向け分子軌道法計算(九大情基センター a , 九州先端研 b ) 本田宏明a , 稲富雄一 a , 真木淳 b

Molecular orbital calculation on many-core processor
using RT parallelization method(RIIT, Kyushu Univ.a , ISITb ) Hiroaki Hondaa , Yuichi Inadomia , Jun Makib【はじめに】
現在、ノートブックからサーバコンピュータに至る多くのコンピュータで複数個の計算コアを持
つマルチコアプロセッサが使用されている。現状では 16 コア以下が大半であるが、今後もコア数は増加してゆき数年以内には非常に多くのコア (> 100 コア) を持つメニーコアプロセッサも利用可
能となると予想されている。また高性能低消費電力化の要請から、今後使用されるハイエンドプロ
セッサは、主要なものでも A: 100 コア以上のメニーコア構成、B: プロセッサあたりの主記憶量の低減ならびに主記憶に対するコアあたりの提供メモリバンド幅の低下、C: 16 ウェイ以上の SIMD
演算器、D: 命令処理におけるインオーダー実行、E: バス方式からパケットスイッチネットワーク
通信型へのコア間通信方式の変更、といった特徴を持つと予想される。
そこで本研究では 2 電子フォック行列 (G 行列) 計算を対象とし、上記のうち A、B 、C の特徴を取り上げ、A と B については RT 並列化法の適用を、C についてはデータ並列計算手法の適用を試み、その有効性の検証を行なう事を目的とした。
【RT 並列化法】
高島らによって開発された G 行列計算に対する㻸㼛㼛㼜㻌㻵㻌㻌㻩㻌㻝㻘㻌㻺RT 並列化法 1 は、元々は各計算プロセッシング 㻸㼛㼛㼜㻌㻷㻌㻩㻌㻝㻘㻌㻵㻰㻵㻶㻔㻝㻨㻩㻶㻨㻩㻵㻕㻘 㻰㻷㻶㻔㻝㻨㻩㻶㻨㻩㻵㻕㻘 㻰㻵㻸㻔㻝㻨㻩㻸㻨㻩㻵㻕㻘 㻰㻵㻸㻔㻝㻨㻩㻸㻨㻩㻵㻕䝕䞊䝍㌿㏦
エレメントに 2.5MB 程度の非常に小さなメモリ㻸㼛㼛㼜 㻶 㻩 㻝㻘 㻵のみ搭載した、分散メモリアーキテクチャのアク㻸㼛㼛㼜 㻸㻩 㻝㻘 㻵㻔㻵㻶㻘㻷㻸㻕 ✚ศィ⟬
セラレーター向けに開発され、図 1 に示すように
㼕㼒 㻔 㻶㻨㻩 㻵 㻒㻒 㻷㻨㻩㻵 㻒㻒 㻸㻨㻩㻵 㻕 㼠㼔㼑㼚
添字駆動型と積分駆動型の G 行列計算方法の混㻳㻵㻶 㻗㻩㻰㻷㻸 㻖 㻔㻵㻶㻘㻷㻸㻕㼕㼒 㻔 㻸㻨㻩㻷 㻒㻒 㻷 㻨 㻵 㻕 㼠㼔㼑㼚
ざりあった計算アルゴリズムとなっている。その㼕㼒 㻔 㻶 㻨 㻵 㻕㻳㻷㻸 㻗㻩㻰㻵㻶 㻖 㻔㻵㻶㻘㻷㻸㻕 㻖 㻞㼑㼘㼟㼑㻳㻷㻸 㻗㻩㻰㻵㻶 㻖 㻔㻵㻶㻘㻷㻸㻕
利点として、各並列計算において必要な G 行列㼑㼚㼐㼕㼒や密度行列 (D 行列) 量は各 4 列ずつとなってお㻳㻵㻸 㻙㻩㻰㻷㻶 㻖 㻔㻵㻶㻘㻷㻸㻕㼕㼒 㻔 㻷 㻨㻩 㻶 㻒㻒 㻶 㻨 㻵 㻕㻳㻶㻷 㻙㻩㻰㻵㻸 㻖 㻔㻵㻶㻘㻷㻸㻕り、データ転送量に関わるコストが著しく少ない㼕㼒 㻔 㻶 㻨㻩 㻷 㻒㻒 㻷 㻨 㻵 㻒㻒 㻸 㻨 㻵 㻕 㻳㻷㻶 㻙㻩㻰㻵㻸 㻖 㻔㻵㻶㻘㻷㻸㻕㼑㼚㼐㼕㼒
点が挙げられる。また、各並列計算部分ではこの㻱㼚㼐 㻸㼛㼛㼜 㻶㻘 㻸
転送された小さなサイズの G,D 行列片による計㻳㻵㻶㻔㻝㻨㻩㻶㻨㻩㻵㻕㻘 㻳㻷㻶㻔㻝㻨㻩㻶㻨㻩㻵㻕㻘 㻳㻵㻸㻔㻝㻨㻩㻸㻨㻩㻵㻕㻘 㻳㻵㻸㻔㻝㻨㻩㻸㻨㻩㻵㻕䝕䞊䝍㌿㏦算が可能であるため、5000 次元の計算を行う場 㻱㼚㼐 㻸㼛㼛㼜 㻵㻘 㻷
合でも 0.6MB 程度のメモリ量にて計算が可能と
図 1: RT 並列化アルゴリズム
なる。
更には 1000 次元程度の場合には 120KB 程度となるため、キャッシュメモリのみにて十分な計算が可能と考えられる。その一方、通常の積分駆動型計算に比較し、2 倍程度積分を重複して計算してしまう欠点を持つ。その結果、コア数の少ない場合では積分駆動型並列計算に比較し実行時間の面で不利だが、数百コア以上の並列計算を行う場合では有利となると予想されている。これらの特性により、数個のメニーコアプロセッサによる並列計算を行う場合では、メニーコアの特徴である
多数のコアの有効利用を可能とし、主記憶バンド幅の低下問題の緩和に役立つ可能性がある。

【データ並列化法による SIMD 計算】
一般に、コンパイラの自動並列化機能を利用してコード中の多重ループ箇所の SIMD 化オブジェ
クトコード生成を行うには、最内ループが容易にループアンローリング可能であることが必要である。しかしながら通常の G 行列計算コードでは、最内ループが縮約計算であるためにそのループ長が外側ループ添字に依存しており、自動的には SIMD 化されない。そこで、ポテンシャルエネル
ギー面作成など、分子座標が異なる複数の入力データを並列に計算する様なデータ並列計算を対象とし、最内ループをデータ並列計算に割り当てることにより SIMD 化オブジェクトコード生成を試みた。また、SIMD 計算を効率的に実行するためには SIMD 演算対象のデータが連続領域に確保されている必要がある。そのため、G 行列計算部分では、積分計算に必要な全ての浮動小数点データ
や D, G 行列の全要素をデータ並列分だけ連続領域に多重化して保持するよう変更を行なった。
このデータ並列計算は前節の RT 並列化法とは全く異なる最適化方法であり、両方の実装を組合せる事が可能である。
【性能評価テスト】
上記の RT 並列化手法ならびにデータ並列化について (ss, ss) 積分を対象とした 逐次実行 G 行
列計算の実装を行い、有効性の検証を行なった。RT 並列化手法において逐次実行する際にはデータ転送部分はメモリコピーにて代用した。使用コンピュータは Intel Xeon X5650 (2.67GHz, 6 コア,
倍精度浮動小数点 2 ウェイ SIMD 演算器) プロセッサ、DDR3 DIMM 16GB メモリからなり、コンパイル条件として Intel コンパイラ 12.1.0 による O3 最適化を使用した。性能計測には Perfmon2
ライブラリ 2 によるハードウエアカウンターデータの取得を行った。
複数の入力ファイルを対象としたデータ並列計算における SIMD 計算の有効利用に対する結果を図 2 に示す。6-31G 基底関数の (H2 O)8 分子を対象とし、分子構造パラメータのみ異なる 32 個の入力データを準備し、1∼32 個の入力に対するデータ並列計算を行なった。グラフ中、折れ線は全
実行浮動小数点命令数中の SIMD 浮動小数点命令数比を示しており、棒グラフはデータ並列を行わないプログラムを基準とした場合の、実行クロックサイクル数/データ並列度 に対する速度向上比を示している。この結果、データ並列が無い場合では、ほぼ 0% であった SIMD 化浮動小数点演算
比が、並列度の増加に従い急増し、16 並列では最大 0.99 と、ほとんどの浮動小数点演算が SIMD化される事が分った。実行クロックサイクル数を計測したデータ並列度あたりの速度向上比はデータ並列度が向上するに従い 2 に近づいており、整数演算等からなる逐次実行の影響もそれほど大きくはなく、全体として 2 ウェイ SIMD 演算器を使用した効率的な計算が実行されていることが分った。この結果から将来のウェイ数の増加に対しさらなる性能向上が可能であると期待される。
RT 並列化についての結果等については当日報告する。

図 2: データ並列計算におけるデータ並列度あたりの速度向上率と全浮動小数点演算数における SIMD 演算数の比
1

H.Takashima et al., J.Comput.Chem., Vol.23, pp.1337-1346, 2002.

2

http://perfmon2.sourceforge.net/

1P129Ab initio 多配置／多参照計算による酸素発生触媒系マンガン複核錯体の電子構造	
 （北大院総合化学 1、北大院理 2）	
  関	 奈々美 1、中山	
 哲 2、武次	 徹也 2	
 Ab initio CASSCF/CASPT2 Calculations for Oxygen Evolving Manganese-oxo Binuclear Complexes（Graduate School of Chemical Sciences and Engineering, Hokkaido Univ.1,
Graduate School of Science, Hokkaido Univ.2）	
 Nanami Seki1, Akira Nakayama2 and Tetsuya Taketsugu2	
 【研究背景】
人工光合成は、現代のエネルギー問題に対する解決策の 1 つとして広く注目を集めており、重要な研究課題となっている。特に近年では、光化学系 Ⅱ	
 (PSII) における酸素発生中心 (OEC)とよばれるマンガン四核錯体が高効率に水分子を分解し酸素分子を生成することから 1、実験・理論の両面から多くの研究がなされ、最近では OEC の詳細な結晶構造が解明されたところである 2。
また、Limburg らは、OEC の部分構造である Mn-(µ-O)2-Mn ユニットを持つマンガン複核錯体：
Mn2O2(tpy)2(H2O)2：tpy = 2,2’:6’,2”-terpyridine が水中で酸素発生を触媒することを報告したが 3, 4、
詳細な反応メカニズムについては未だ解明されていない。
遷移金属元素や大きな配位子を含むような、規模の大きな系に対して理論計算を行う場
合、多くは DFT 法が用いられる。DFT 法は、適度な計算コストで分子の構造や励起エネルギーなどを実験とも比較可能な精度で得られる非常に有用な計算手法であるが、単配置理論であるため、開殻系の非制限計算においてはスピン混入がおこり、正しい電子状態を解くことができない。
一方で、ab initio 法である MCSCF 法は、波動関数が電子配置の線形結合で表される多配置理論であり、スピンの固有関数になる正しい電子状態を解くことができるが、遷移金属が複数含まれる系では活性空間が増大して計算コストが高くなるため、規模の大きな系にはほとんど適用されてこなかった。しかし、酸化還元反応を記述するにあたっては、各反応中間体における中心金属の酸化数、すなわち開殻 d 電子の数に応じたスピン状態の変化を正確に追うことが不可欠であることから、多配置理論に基づく電子状態計算を行うことが求められている。
以上の背景をふまえ、本研究では、Limburg らが報告したマンガン複核錯体に対して多配置／多参照計算を適用し、多配置波動関数に基づく電子状態を求めることを試みた。中心金属のスピン状態に依存する物理量の一つである有効交換積分 (J) を DFT 法・CASSCF 法を用いて算出し、酸化数に従った変化と計算手法による違いについて議論する。

【計算方法】
マンガン複核錯体：Mn2O2(tpy)2(H2O)2 を計算対象とした。はじめに、2 つのマンガン原子の酸化数の組み合わせを III-III、IV-III、IV-IV、V-IV の 4 通り設定し、各系の最高スピン状態について DFT(UB3LYP) 法を用いて構造最適化を行った。ここで、分子構造について等原子価状態では C2h 対称性、混合原子価状態では Cs 対称性の拘束を用いた。得られた分子座標を用いて、DFT(UB3LYP) 法で最低スピン状態の 1 点計算を行い（1 重項の場合は broken symmetry 法を適用）、さらに CASSCF/CASPT2 計算を行った。CASSCF 計算における活性空間には、2 つのマンガンの d
軌道を主成分とする 10 個の軌道を含めた。CASPT2 計算においては、コストを削減するため下から 100 個の軌道を frozen core として扱い、60 個の軌道からの励起を考慮した。基底関数は、マンガンに対し Stuttgart/Köln 擬ポテンシャル、その他の原子に対し Sapporo-DZP-2012 を適用した 5。
DFT 計算は Gaussian 09、CASSCF/CASPT2 計算は MOLPRO 2008, 2010 を利用して行った。

【結果と考察】各 酸 化 数 の 組 み 合 わ せ で 最 適 化 構 造 を 比 較 す る と ( 表 1) 、 マ ン ガ ン 原 子 間 距 離(d[Mn–Mn]) は酸化数の増加に従って大きくなり、配位結合距離 (d[Mn–N]、d[Mn–O(H2O)]) は小さくなる傾向が示された。ただし、V–IV 状態については、価数による差がほとんど見られなくなった。また IV–IV 状態に対する最高スピン・最低スピン状態の相対エネルギーを比較すると(表2)、ともに低スピン状態が安定で J が反強磁性相互作用を示す点は一致したが、DFT と CASSCF
で J の絶対値に 102 のオーダーの差が生じた。今回の計算から得られた CASSCF 自然軌道には架橋酸素の p 軌道の成分が混合してきていることから、活性空間における架橋酸素を通じた電子の相互作用は十分表現されており、J 値の差は波動関数の違いによって生じたものと考えている。
詳細は当日報告する。

表 1. 主な原子間結合距離（Å）
酸化数

d[Mn–Mn]d[Mn–N]:longd[Mn – N]:shortd[Mn – O(H2O)]

III – III

2.740

2.331

2.173

2.169

IV – III

2.746(IV)2.062 / (III)2.273(IV)2.028 / (III)2.079(IV)2.159 / (III)2.195

IV – IV

2.752

2.056

1.995

2.124

V – IV

2.775(V)2.077 / (IV)2.077(V)2.013 / (IV)2.013(V)2.118 / (IV)2.118

表 2. IV–IV 状態に対する相対エネルギー（kcal mol-1）
SeptetSinglet

DFT

2.539

0

J = -74.01 cm-1

CASSCF

0.026

0J = -0.77 cm-1

【参考文献】
1. V. K. Yachandra, K. Sauer, M. P. Klein, Chem. Rev. 96, 2927 (1996).
2. Y. Umena, K. Kawakami, J.-R. Shen, N. Kamiya, Nature 473, 55 (2011).
3. J. Limburg, J. S. Vrettos, L. M. Liable-Sands, A. L. Rheingold, R. H. Crabtree, G. W. Brudvig, Science
283, 1524 (1999).
4. J. Limburg, J.S. Vrettos, H. Chen, J.C. de Paula, R.H. Crabtree, G.W. Brudvig, J. Am. Chem. Soc. 123,423 (2001).
5. T. Noro, M. Sekiya, and T. Koga, Theo. Chem. Acc. 131, 1124 (2012).

1P130摂動法に基づいた電流存在下での電子構造計算手法の研究（京都大院・工）池田裕治，瀬波大土，立花明知

Study of Perturbative Approach for Calculating Electronic Structures under
the Existence of Electric Currents(Kyoto Univ.) Yuji Ikeda, Masato Senami, Akitomo Tachibana【序】
単分子トランジスタに代表されるようなナノ材料の電気伝導特性の理論的解析としては，非平衡 Green 関数法のような手法による，デバイス全体としての巨視的な電流電圧特性の計算などが主流である．しかし，ナノ材料の電気伝導特性を本質的に理解し，材料設計などに活かすためには，巨視的な電気伝導特性のみでは不十分であり，系の波動関数を基に，局所的な電気伝導特性の解析を行わなければならない．そのためには，電気伝導状態の波動関数を求める手法が必要となる．そこで本発表では，電子間の交換相互作用を取り入れた，電気伝導状態における摂動波動関数を構成する手法について報告する．また，Rigged QED において定義される局所電気伝導率テンソルを求め，ナノ材料の局所的な電流電圧特性についての解析結果について報告する．
【計算方法及びモデル】
左右電極によるチャネル部への効果は電子のエネルギー E に依存する自己エネルギー行列 ΣL,R (E)により取り入れることができる．そのため，ΣL,R (E) をチャネル部の状態密度 w(E) と掛け合わせてエネルギー領域で積分することにより，電気伝導状態に対応する摂動ハミルトニアンを得る
ことができる．すなわち，∫ ∞(1)
Σ(E)w(E)dE,
H =−∞である．ここで，Σ(E) は ΣL,R (E) からなり，

Σ(E) = (ΣL (E) + Σ†R (E))(1 − fL (E))fR (E)+ (Σ†L (E) + ΣR (E))fL (E)(1 − fR (E)),
である．fL,R (E) はそれぞれ左右電極の Fermi-Dirac 分布関数である．ΣL,R (E) は 1 次元電極に対する Green 関数を用いて計算される [1]．ΣL,R (E) はチャネル部からの電子の流出に対応する．一†方，チャネル部への電子の流入は，計算の簡略化のため ΣL,R (E) で表現されるとしている．w(E)は，w(E) =γp1∑,
π p (E − ϵp )2 + γp2と表される．ここで，ϵp は軌道 ϕp のエネルギーであり，γp は ϕp の準位の広がりである．ϵp , γp は摂動後の軌道から決定されるため，H (1) に依存している．そのため，ϵp , γp と H (1) とが互いに無矛盾となるまで摂動計算を繰り返す．
H (1) による波動関数の摂動は，電子間の交換相互作用を取り入れるため，Coupled Perturbed
Hartree-Fock (CPHF) 方程式によって求める．H (1) は一般に非エルミート行列であるため，以下のように取り扱い，対応する波動関数を求める．まず，H (1) = H H(1) + H A(1) ,とし，H (1) をエルミート行列 H H(1) と歪エルミート行列 H A(1) とに分ける．そのそれぞれに対し，摂動後の軌道 ϕp は，変換行列 U を用いて，
∑
ϕp =
ϕ(0)
q Uqp ,qと書ける．U は生成子 κ により U = exp(κ) と表され，κ は CPHF 方程式によって求められる．
ここで，H H(1) に対しては κ が歪エルミートとなるのに対し，H A(1) に対しては κ はエルミート

ˆとなる．得られた波動関数を用いて，局所電流密度ベクトル ⃗j(⃗r) が計算される．また，Rigged
QED において，ˆ⃗ˆ r),⃗ˆj(⃗r) = ↔
σ ext (⃗r)D(⃗

y
x

↔
ˆ

H

z
のように定義される局所電気伝導率テンソル σ ext (⃗r)[2] が計算される．
S
本研究では主にベンゼンジチオール及びその誘導体に対しての計算を行う．図 1 に，本研究で使用するベンゼンジチオールの構造を示す．
構造は Hartree-Fock レベルで最適化されている．基底関数は STO-6G
Cを用いる．無摂動状態における電子構造については Hartree-Fock 法
により求める．電流存在下での電子構造は独自に開発したプログラム
を用いて計算される．両端の S 原子は，3px , 3py 軌道を介してそれぞ 図 1: ベンゼンジチオール
れ −2.0 eV の結合定数を持って 1 次元 Au ナノワイヤ電極に接続されているとする．
【結果と考察】本抄訳では，ベンゼンジチオールの局所電流密度ベクトルの解析結果を示す．図 2 に V = 1.50 Vのバイアス電圧印加時における局所電流密度ベクトル ⃗j(⃗r) を示す．電流密度の大きさや向きが空間各点で大きく異なっており，電気伝導特性の空間依存性が存在していることが確認できる．また，y = 0 bohr 面よりも y = 1.651 bohr 面の方がより広い領域で電流密度が大きく，電流が π 軌
道を通って流れていることが確認できる．図 3 に V = 2.00 V のバイアス電圧印加時における局所
電流密度ベクトル ⃗j(⃗r) を示す．V = 1.50 V のバイアス電圧印加時と比較して，全体的に電流密度が大きくなっていることが確認できる．今回の計算手法では V = 1.00 V 以下では電流密度は生じなかった．
本講演では，ベンゼンジチオールの誘導体に対する結果との比較，及び局所電気伝導率テンソ
ルについても発表する．
【参考文献】[1] T. Shimazaki, H. Maruyama, Y. Asai, and K. Yamashita, J. Chem. Phys. 123, 164111 (2005).
[2] A. Tachibana, J. Mol. Model. 11, 301 (2005); J. Mol. Struct.: THEOCHEM 943, 138 (2010).
−3

8

−3

1.0 × 10 8
SC0

z (bohr)HCz (bohr)z (bohr)C−40

1.0 × 10
C
H

C

−4
0

1.0 × 10−4

1.0 × 10

H

C

C

S

S

−8
−8

1.0 × 10
S

C

H

−3

1.0 × 10 8

0

8

−5
1.0 × 10−8−80

x (bohr)x (bohr)

(a)

(b)

8

−5
1.0 × 10−8−8

1.0 × 10−508y (bohr)(c)

図 2: V = 1.50 V のバイアス印加時の局所電流密度ベクトル ⃗j(⃗r) (a.u.)．(a) y = 0 bohr 面．
(b) y = 1.651 bohr 面．(c) x = 0 bohr 面．
8
1.0 × 10−388
1.0 × 10−3

S

C

C

C

C

z (bohr)HCz (bohr)z (bohr)

S

C

H

0
1.0 × 10−4

0
H

C

0
1.0 × 10−4

H

S

0

1.0 × 10−4
C
S

−8
−8

1.0 × 10−38−5
1.0 × 10−8−80

x (bohr)x (bohr)

(a)

(b)

8

−5
1.0 × 10−8−8

1.0 × 10−508y (bohr)(c)

図 3: V = 2.00 V のバイアス印加時の局所電流密度ベクトル ⃗j(⃗r) (a.u.)．(a) y = 0 bohr 面．
(b) y = 1.651 bohr 面．(c) x = 0 bohr 面．

1P-131
酢酸-メタノール錯体における二重プロトン移動反応についての理論的研究（名工大院・工） 藤澤佑紀、志田典弘

Theoretical study of double proton transfer reaction of acetic acid -methanol complex(Nagoya Institute of Technology) Yuki Fujisawa, Norihiro Shida【序】 プロトン移動反応は化学反応の代表的な反応の一つであり、特に多重プロトン移動反応は、生体内における諸反応で重要な役割を果たしている。しかし、その反応メカニズムについては未だ不明な部分が多く、理論と実験の両面から様々な分子を用いた研究が行われてきた。多重プロトン移動反応のモデル系としては、ギ酸二量体（FAD）が広く用いられてきたが、今回我々は、FAD よりもより対称性の低い酢酸-メタノール
錯体について、量子力学的計算による理論的側面から、より一般的な反応メカニズムについての解析を行った。プロトン移動反応では、プロトンの軽さから生ずる量子効果が重要な役割を演じる。この効果には、反応障壁を通り抜けるトンネル効果や、複数の経路の存在を示す多次元効果が含まれる。トンネル効果については一般にトンネル二準位分裂（TS）がその大きさを知るための指標となる。しかし、この酢酸-メタノール
錯体の反応障壁（16.14 kcal/mol）は、プロトン移動反応の反応障壁としてはかなり大きく、TS は実質的には 0 に等しい。そこで本研究では、反応曲面法を用いて時間依存の波動関数の振る舞いを直接調べること
により、反応メカニズムの解析を行った。
【計算手法】 多次元の反応曲面を定義するには、q1~q4 の４つの対称座標を反応座標に用いた。他の自由度については、逐次エネルギーが最低になるよう最適化し、反応曲面を定義した。
(1.a)
(1.b)
(1.c)
図１：内部座標(1.d)

q1~q4 の座標はそれぞれ、「２つのプロトンの同期的な運動」、「酢酸分子とメタノール分子間の対称な伸縮運動」、「２つのプロトンの非同期的な運動」、「酢酸分子とメタノール分子間の非対称な伸縮運動」を表す。
なお、q1, q4 に関しては、以下の対称性を満足するように座標変換を施したものを用いた。
(2)電子状態の計算には HF 法と CCSD 法を用い、まず HF/6-31G**による制限付き構造最適化を独立な９６０点で行うことで反応曲面を定めた。次に、その各点について CCSD/6-31G**による計算により、ポテン
シャルエネルギー曲面（PES）を求めた。各点におけるエネルギーは、次の解析関数により、最小二乗法を用いてフィットした。
(3)
量子力学的計算には、四次元の反応曲面 Hamiltonian の固有解を、MCSCF 法的な基底関数展開法により下から約 10000 個計算した。
(4)()()
時間依存の波動関数は、これらの定常解の疑似的な完全系で展開することにより表現した。
(5)

【PES】 得られた平衡状態と遷移状態の対称座標は、以下の値となった。（単位: a.u.）
平衡状態 (min):
遷移状態 (ts):
図２~図５は、反応曲面上の PES を、様々な部分空間で切り出して等高線図で表したものである。

図 2. PES q1 vs. q2 (min)図 3. PES q1 vs. q3 (min)図 4. PES q1 vs. q3 (ts)

図 5. PES q1 vs. q3 (q2: 最適化)

図２は平衡状態での２つのプロトンの同期的な運動(q1)と分子間の対称な伸縮運動(q2)の自由度に対応した PES であり、対称二極小の形をしている。図３、図４は、平衡状態と遷移状態での、２つのプロトンの同期的な運動(q1)と非同期的な運動(q3)の PES である。q2 の値の変化に伴い曲面の形状も大きく変化し、q1 とq3、q2 の間には大きな相関があることが窺える。図５は、q2 の自由度を全て最適化した q1 対 q3 の曲面であり、
丁度図 3 と図 4 の中間くらいの形状であることが分かる。
【結果】 図６は、PES の左側の極小値（min）に局在化した初期波束の反応の初期段階（t=250 a.u.）における振る舞いを、q1, q3 の関数として表したものである。なおこの計算では、全エネルギー（4257 cm-1）を反応障壁（5624 cm-1）よりも低く設定し、反応過程が 100% トンネル効果となるようにした。

図６. 波束の時間発展 q1 vs. q３

図７. q1=0.0 における波束の振る舞い

図８. q1=3.0 における波束の振る舞い

図６からは、q3=①1.9, ②1.0, ③0.0, ④-1.2 付近を中心とした右側へ伸びる 4 本のピークが確認される。
このうち①, ②, ④のピークは非同期的な運動に対応し、③は同期的な運動に対応する。また相対的には、③, ④に比べ①, ②が大きな確率振幅を持つピークであることがわかる。図７、図８は、与えられた q1 における波束の確率振幅を q3 の関数として時間の経過と共にどのように変化したかを表したものである。ここで図７は反応の分水嶺である q1=0.0 の振る舞いを表したもの、図８はほぼ反応が完了した q1=3.0 の振る舞いを表したものである。図７より、反応の分水嶺では時間の経過に関わらず、常に①, ②の２つのピークが優勢となることがわかる。これに反し図８では、確率振幅は小さいものの、③の同期的な運動のピークの存在が読み取れる。この q1 の値によって新たに出現したり消滅したりするピークの解釈については、現在解析を進めている。しかしながら図７、図８全体を見ると、大きな確率を持つのは常に①, ②付近を中心とした非同期的な運動のピークであることがわかる。講演では、これらの詳細について報告する。

1P132ホモロジーモデル内部の水分子位置決定に関する理論的研究（九大先導研）緒方龍展, 蒲池高志, 吉澤一成

Theoretical Study on Prediction of Water Molecules’
Locations in Homology Models(IMCE Kyushu University) Tatsunobu Ogata, Takashi Kamachi, Kazunari Yoshizawa【序】
水分子はタンパク質のフォールディング、構造安定性、機能などに重要な役割を果たすた
め、水分子の情報はタンパク質の計算モデルに不可欠である。また、タンパク質の構造決定
法には X 線結晶構造解析、中性子線回折、NMR などがあるが、膜タンパク質のような結晶化が難しいタンパク質にはそれらの手法を用いることができないため理論的手法による構造決定が必要である。タンパク質の理論的構造決定法の一つにホモロジーモデリングがある。
この手法では、目的のタンパク質とアミノ酸配列が類似し、かつ、構造が既知であるタンパ
ク質を鋳型として目的のタンパク質のアミノ酸残基を並べ、構造を決定する。鋳型の配列類
似度が 30%以上であれば主鎖レベルで 2 Å の誤差での予測が可能である。しかし、ホモロジ
ーモデリングでは水分子の位置を決定することができない。本研究ではホモロジーモデリン
グにより作成したタンパク質の内部の水分子位置を決定する手法を開発した。
【計算方法】
ホモロジーモデリングには Discovery Studio を用いた。鋳型となるタンパク質の構造情報
を PSI-BLAST により検索し、配列類似度が高いものから 3 つを選択し、鋳型を作成した。鋳型の作成、モデルの構築には MODELER を用いた。アミノ酸残基のプロトン状態決定にはGeneralized Born 法を用い、PROPKA を用いてクロスチェックを行った。なお、2 つの方法で結果が異なる際には各アミノ酸残基の一般的なプロトン状態（HIS(0)、GLY(-1)、ASP(-1)、
ARG(+1)、LYS(+1) ; ( )内は電荷）とした。分子動力学計算には力場として CHARMM22 を用い、計算プログラムには NAMD を用いた。開発した手法を評価するために、水の位置が精度よく決定されている高解像度の X 線結晶タンパク質を対象とした。開発した手法により決定された水分子が X 線結晶構造解析により得られている水分子と対応するかを評価した。
また、水分子を除いた X 線結晶構造とそのホモロジーモデルについてそれぞれ水分子を予測することで、タンパク質構造の精度が与える影響について確認した。
【結果と考察】本研究ではタンパク質内部の空隙に存在する水分子(buried water)の位置を分子動力学シ
ミュレーションにより予測する手法を開発した。buried water はタンパク質の構成原子により外部へ到達することができず、かつ、溶媒露出度が 0 である水分子である。我々はタンパ
ク質について原子間距離を計算することにより空隙を探索し、そこに水分子を挿入するプログラムを開発した。作成した構造について分子動力学シミュレーションを行った（50 ps の平衡化、0.001 kcal/mol·Å の収束条件での共役勾配法による構造最適化）。Li および Jan らの報告によると[1]、buried water はタンパク質との非結合性相互作用により 12 kcal/mol 以上安定化されている。そこでシミュレーションの結果得られた構造から各水分子の安定化エネ
ルギーを計算し、-12 kcal/mol 以上の水分子を削除した。不安定な水分子がなくなるまでシ
ミュレーションと削除を繰り返すことによりタンパク質内部の水分子の位置を決定した。水分子が一つのみ存在する空隙から複数の水分子が存在する空隙まで水分子位置を精度よく決定できた（図１）。得られた結果を表１に示す。X 線結晶解析により得られている水分子に対し、水分子を除いた結晶構造に対する水分子位置の予測精度は非常に高く、また、ホモロジーモデルについては RMSD 値が良ければ水分子位置を高精度で決定できることが判明した。これにより本手法は対象とするタンパク
質の構造精度が高ければ水分子を正しく予測するこ

図１．X 線結晶解析から得られているとが可能であり、ホモロジーモデルに対する水分子

水分子（青球）と本手法で予測した水位置決定手法としても有効であることが示された。

分子（赤球）

表 1．X 線結晶構造解析で決定された buried water の位置と本手法によって決定した水分子の位置の比較結果予測した水の個数
PDB ID

解像度(Å)

buried watera

RMSD(Å)b

X 線結晶

ホモロジーモデル

0.80 Åc

1.00 Åc

0.80 Åc

1.00 Åc

3A5Q

1.80

14

0.64

13

13

8

9

2DE9

1.30

13

0.74

12

12

9

9

3PTL

1.30

7

0.90

6

6

6

7

3KBJ

2.00

12

1.27

12

12

8

8

3O6A

2.00

15

1.63

11

12

11

11

5PTP

1.34

10

1.73

7

7

5

8

1TVN

1.41

9

2.37

9

9

4

4

70

71

51

56

合計

80

a X 線結晶構造解析で得られた buried water の個数
b X 線結晶構造に対するホモロジーモデルの誤差
c X 線結晶構造解析で観測された水分子位置と本手法により決定した水分子位置の誤差【参考文献】[1] Li, Z.; Jan, H. PROTEINS 1966, 24, 433.

1P133
時間依存シュレディンガー方程式に基づく
量子シミュレーションの GPU による高速化(豊橋技術科学大学院）

増田拓弥, 関野秀男

High Performance Simulation of TDSE using GPU(TOYOHASHI UNIVERSITY OF TECHNOLOGY)
Takuya Masuda， Sekino Hideo【研究背景・目的】
原子や分子の動的物性の予測を行う場合，量子力学における基本方程式であるシュレディンガー方
程式を解くことは必須である．しかしながら，分子の電子状態に対する予測の精密な結果を得るためには膨大な計算領域と解析時間を有するため，高速演算のアルゴリズム開発と並列環境の実装は重要である．
並列演算に特化した演算装置として GPU(GraphicsProcessingUnit)がある．GPU は，主に画像処理用の演算装置であるが，近年，その高い演算処理能力から， 画像以外の一般的な用途に利用する
GPGPU(GeneralPurposeGPU)という概念が注目されている．
そこで，本研究は NVIDIA 社製 GPU と C 言語の統合開発環境 CUDA(Compute Unified device
Architecture)を用いることによる時間依存シュレディンガー方程式に基づく量子シミュレーションの高速化を行った．

【シュレディンガー方程式の差分化】
量子力学の基本方程式として，2 次元のシュレディンガー方程式があり，次のように記述される．

2  22i  x, y, t   
 2  x, y, t   2  x, y, t 式(1)t2m  xy
また，流体の波動現象を記述するものとして，次の波動方程式がある．
2

2
2
2 




z
x
,
y
,
t

c
z
x
,
y
,
t

z x, y, t  式(2)

2
2
2
t
y
 x

式(2)にテイラー展開を適用して解くと，現在と過去の変位量を用いて未来の変位量を推定する次の方程式を導くことができる．

z x, y, t  1 

4  8c 2t 2 d 2
t  22c 2t 2 d 2zx  1, y, t   zx  1, y, t   zx, y  1, t   zx, y  1, t 式(3)z  x, y , t  z x, y, t  1 t  2t  2t  2式(1)(2)の比較から，シュレディンガー方程式は流体の波動方程式との類似性が高いことがわかる．
つまり，流体理論で適用できることは量子理論でも同様に適用できることが考えられるため，式(1)を式(3)の形式に解くことで，量子の動的物性の計算における GPU による高速化の影響を比較できる．
量子の動的物性のモデル化に際し，2 次元断面は n×n のグリッド上に格子点を横間隔  x，縦間隔 y で配置されたメッシュで表現し，各格子点における波動関数の値は式(1)を用いて推定する．
ψ(x,y,t)は量子の波動関数であり，x,y は x，y 方向の離散化された空間座標，t は時間間隔 t の離散化された時間座標である．導関数の近似にテイラー展開による 2 階差分近似を用いると，x、y に関する 2階導関数を得ることができ，中心差分近似を用いると t に関する 1 階導関数を得ることができる．これ
らの導関数を式(1)の微分項に代入し，整理すると次のように書くことができる．
i  x  x, y, t    x  x, y, t   2 x, y, t   x, y  y, t    x, y  y, t   2 x, y, t  x, y, t  t    x, y, t  t   t 式(4)
m
x 2
y 2

この方程式は，前回・現在のその格子点における波動関数値と，現在の隣接格子点における波動関数値をもとに，時間 t だけ経過した未来の波動関数値を推定するものである．

【CUDA を用いたプログラミング】
GPU は，数万スレッドにも及ぶ CPU より圧倒的に多い数のマルチ・スレッド方式で処理することに

よって，高い並列演算処理を実現することができる．
この GPU を動作させるために必要なものが，C 言語の統合開発環境 CUDA である．プログラム構成
は，CPU 側で動き，配列の確保などを行うホストコードと，GPU に処理させたい内容を記述するデバ
イスコードで成り立っている．メモリ構造には，高速アクセスが可能なメモリであるシェアードメモ
リと，大容量であるがアクセスが低速なグローバルメモリがある．また，グリッド，ブロックという階層を用いることで，膨大なスレッド数を管理している．つまり，GPU の性能を発揮させるためには，使用するメモリの特性とスレッド階層構造を考慮したプログラミングが求められる．
本研究では，GPU を使用したプログラムとして，グローバルメモリのみを使用したプログラム，グローバルメモリとシェアードメモリを併用したプログラムの 2 つを作成した．
グリッド内ブロック数，
ブロック内スレッド数はともに 2 次元配置とした．1 格子点当たり 1 スレッドで計算するようにし，その中で時間間隔 t 毎にデバイスコードに記述した式(3)の処理を時間発展させることで，時間依存シュ
レディンガー方程式を解いた．

【検証方法】
CPU，ブロック内のスレッド数を 16×16 で統一した GPU(グローバルメモリ使用)，GPU(シェアー
ドメモリ使用)によってシュレディンガー方程式を計算する 3 つのプログラムを用いて，スレッド数
128×128，256×256，512×512，1024×1024，2048×2048，4096×4096 における各プログラムの計算
時間を計測する．なお計測は，それぞれ 3 回ずつ行い，その平均値をとった．

【結果】
スレッド数毎の各プログラムにおける計算時間の計測結果を表 1，CPU を 1 とおいた場合の計算時間比を図 1 に示す． なお，図 1 の縦軸は計算時間比であり，[CPU 計測時間/GPU 計測時間]とした．

表 1

各プログラムにおける計算時間(単位：ms)

128*128

256*256

512*5121024*1024

1206.705014.2020691.59

85877.56

GPU(global memory)

50.26

147.15509.78

1987.277478.6928828.77

GPU(shared memory)

31.88

88.01

345.381357.174913.50

18697.55

CPU

図 12048*2048

4096*4096

347234.73 1585193.70

CPU に対する GPU の計算時間比【考察】
表 1 より，計算速度は GPU(シェアードメモリ)>>GPU(グローバルメモリ)>>CPU となった．つまり，
シュレディンガー方程式の計算の高速化において，GPU はともに CPU より有効で，メモリとしてシェ
アードメモリを使用することで最も高速化されることがわかった．
また，図 1 からは，スレッド数の拡大に伴って，計算の高速化がなされていることがわかる．つまり，GPU の性能は，より膨大なスレッド数での計算において最大限に活かされる．
シュレディンガー方程式と流体方程式の GPU による計算結果の比較も行い，複素数が高速化にどのような影響を与えているか考察する．

【参考文献】[1] 青木尊之, 額田彰. はじめての CUDA プログラミング. 工学社, 2009

1P134
外場中に置かれた原子系の Rigged QED シミュレーション(京大院・工) 宮里敏秀，瀬波大土，池田裕治，高田崇二郎，立花明知

Rigged QED Simulation of Atomic Systems in an External Field(Dept. of Micro Eng., Kyoto Univ.)Toshihide Miyazato, Masato Senami, Yuji Ikeda, Sojiro Takada, Akitomo Tachibana

序
量子的な現象の解明に多くの成果を挙げてきた量子電磁力学 (QED) であるが、そこでは電子と光子のみが量子化され、原子核を含む束縛系の扱いには十分でない。一方で、近年半導体分野に於いて注目されている high-k 材料などでは原子核が大きな寄与を持つことが分かっており、原子核をも統一的に扱える理論に基づく研究が望まれている。そこで本研究室では、これまで立花の提唱する Rigged QED[1] に基づいて、静的な局所誘電物性について研究を行って来た [2]。
本研究では、Rigged QED に基づき、量子化された電子場、原子核の場、光子場の時間発展をシミュレートするために開発した計算コード・QEDynamics[3] を用いて、電子場、原子核場と光子場との相互作用を時々刻々計算することで、それらの時間発展の様子を視覚化し、動的な局所誘電物性について議論を行う。今回は、特に対象とする系を取り囲む環境の作り出すスカラーポ
テンシャルに対する系の応答についてシミュレーションを行った。
理論と計算方法
本研究では、輻射光子場は無いものとして相互作用光子場のみを扱う。また、ベクトルポテン
シャルの影響は微小であるものとして無視する。以下ではこれらの仮定の下での式を示す。
電子場の演算子 Ψ̂、原子核 a の場の演算子 χ̂a 、光子場の演算子 Â0 を、電子の生成・消滅演算子 ên 、原子核の生成・消滅演算子 fˆan を用いて次のように表現する。
∑∑
Ψ̂(t, ⃗r) =ψn (⃗r)ên (t), χ̂a (t, ⃗r) =χan (⃗r)fˆan (t), Â0 (t, ⃗r) = Â0A (t, ⃗r) + Â0M (t, ⃗r).
nn

但し系の内部でやり取りされる光子場を表す Â0A は電子・原子核の生成・消滅演算子により、系と周囲環境の間でやり取りされる光子場を表す Â0M は環境の電荷密度演算子 ρ̂ext により展開される。(Za :原子核 a の電荷)∫∫†s)Ψ̂(t, ⃗s) ∑χ̂†a (t, ⃗s)χ̂a (t, ⃗s)
3 Ψ̂ (t, ⃗Â0A (t, ⃗r) = − d ⃗s+
Za d3⃗s,|⃗r − ⃗s||⃗r − ⃗s|a∫ρ̂ext (t, ⃗s)Â0M (t, ⃗r) = d3⃗s.
|⃗r − ⃗s|
特に本研究では、環境のケットベクトル |M ⟩ に対しρ̂ext (t, ⃗r) = −σ̂(t)δ(z − d) + σ̂(t)δ(z + d),σ̂(t)|M ⟩ = σ(t)|M ⟩.

とする。これは、系が空間的に一様な z 軸方向を向いた外部電場中に置かれていることに相当する。
電子場、原子核場の演算子は以下の方程式を満たす。(ma :原子核 a の質量)
]
[
]
[
∂
1
∂
1
2
2
(−i∂i ) + Za Â0 χ̂a .
i Ψ̂ =(−i∂i ) − Â0 Ψ̂, i χ̂a =
∂t
2
∂t
2ma
尚本研究では、環境の電荷密度分布は系の状態に影響されないとしている。
この方程式を解くことにより、電子場・原子核場の時間発展を追うことが出来る。本研究では
特に、電子密度と、環境との相互作用によって変化した電子分布によって作られる局所分極ベク
トルの振舞いについて議論する。
∫
1⟨Ψ̂† (t, ⃗s)Ψ̂(t, ⃗s)⟩ˆ⟨ρˆe ⟩ = ⟨Ψ̂† Ψ̂⟩, ⟨P⃗e ⟩ = − grad d3⃗s.
4π|⃗r − ⃗s|

結果と考察環境の電荷密度分布が時間的に一定である場合の、水素原子についてのシミュレーション結果を示す。図 1 は、電子密度の初期状態からの変化を表す。図 1(b) を見ると、静電場中に置かれた
場合でも電子密度が振動していることが分かる。その振動周期は、遷移の起こっている準位間のエネルギー差に相当している。
また、環境との相互作用によって変化した電子分布によって作られる局所分極ベクトルの分布の様子を図 2 に示す。従来の巨視的な誘電物性論では得られない、原子内部での局所的な分極の様子が可視化できた。
5.0E-073.0E-06
1.5E-06

0.0E+00

0.0E+00-5.0E-07-1.5E-06-1.0E-06-3.0E-06-1.5E-06-2-1
y (bohr)

20

1-2.0E-06

0

1
2

-1

z (bohr)-2(a) 時刻 t=2.0[a.u.] における yz 平面での変化-2.5E-06
0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00
time[a.u.](b) 点 (x,y,z)=(0.0,0.0,0.8)[a.u.] における時
間変化

図 1: 電子密度の初期状態からの変化

2.0
4.0E-07

1.5

z [bohr]

1.0

2.0E-070.5
0.0

0.0E+00-0.5-2.0E-07

-1.0
-1.5

-4.0E-07

図 2: 環境との相互作用により変化した電子によって作られた、局所分極ベクトルの yz
平面上での分布(カラーバーが電子の作り出す静電ポテンシャルを表し、その勾配である図中の矢印が局所分極ベクトルのyz 平面への射影を表す。)-2.0-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0
y [bohr]以上の結果に加えて発表では、環境の電荷密度分布が時間に対し振動的であった場合、及び対象となる系が 2 原子系の場合についても報告する。
結論
Rigged QED によって、場の時間発展の様子を時々刻々追えることが確認された。また局所的に定義された物理量を計算することで、従来の巨視的な誘電物性論では見ることのできなかった、
より局所的で詳細な応答が見られることが分かった。今後は対象とする系や環境を拡大し、より様々な状況における動的で局所的な誘電物性について調べる予定である。
参考
1 A. Tachibana, J. Mol. Mol. THEOCHEM. 943, 138 (2010); A. Tachibana, J. Mol. Modelling 11, 301 (2005)
2

K. Doi et al., Jpn. J. Appl. Phys. 47, 205 (2008); A. Fukushima et al., Jpn. J. Appl. Phys. 49, 121504 (2010);
A. Fukushima et al., Jpn. J. Appl. Phys. 49, 111504 (2010) ; M. Senami et al., Jpn. J. Appl. Phys. 51, 031101(2012)
3
QEDynamics, M. Senami, K. Ichikawa, A. Tachibana (http://www.tachibana.kues.kyoto-u.ac.jp/qed)

1P-135
量子化学計算を用いたラジウムの水和構造探索(お茶大・人間文化創成科学 1，産総研 2，立教大 3)	
 松田 彩 1，古明地 勇人 2，望月 祐志 3，森 寛敏 1

Theoretical investigation of the hydration structure of Ra(II)using quantum chemical calculations(Ochanomizu Univ. 1, AIST2, Rikkyo Univ. 3)
Aya Matsuda1, Yuto Komeiji2, Yuji Motizuki3, Hirotoshi Mori1【序論】
	  金属イオンの水和構造は，金属イオンの基礎的なプロパティの一つであり化学的に重要な情報である．アルカリ土類金属の最重元素であるラジウムは，その発見から100 年経っているにも関わらず，実験的・理論的な研究の数は限られ，水和構造に関するデータが得られていなかった．本研究の目的は，量子化学計算手法による Ra2+ の水和構造の調査である．Ra2+ の水和構造を，他アルカリ土類金属イオンと比較しつつ考察したので報告する．
【計算方法】
	  以下の静的水和モデルおよび動的水和モデルにより Ra2+ の水和構造を検討した．
１．静的水和モデル
	  第二水和圏以降の溶媒効果を連続誘電体モデル（PCM 法）によって記述した水和
クラスターモデル [Ra(H2O)n]2+（n = 1-9）を構築し，HF・B3LYP・MP2 レベルの電子状態計算よる構造最適化・振動数解析から，水和物の熱力学的安定性の評価を行った．初期構造は MMFF94x 力場を用いた古典分子動力学計算結果より，スナップシ
ョット構造を抽出することで作成した．他のアルカリ土類金属イオン（Mg2+, Ca2+, Sr2+,
Ba2+）の水和についても同様な評価を行い Ra2+ 水和系と比較した．Natural Population
Analysis（NPA）, Natural Bond Orbital（NBO）, Localized Molecular Orbital Energy
Decomposition Analysis（LMO-EDA）[1] によって，アルカリ土類金属の違いによる，
水−金属イオン間結合特性の差異を検討した．
２．動的水和モデル
	  Ra2+ を中心にした半径 9.0Å の水球モデルを構築し，FMO-MD [2] 法により，動径分布関数および角度分布関数を解析した．FMO 計算のレベルには，3 体効果までを取り扱った FMO3-HF を用い，温度 300K，時間刻み 1 fs/1 step の条件で MD 計算を実施した．定温アルゴリズムには Nosé-Hoover Chains 法を用い，1 ps に亘る平衡化を行った後，続いて 2ps のプロダクションランを実施し，その結果を統計解析に用いた．
上記計算において，アルカリ土類金属の基底関数には MCP 基底：MCPdzp-CV [3] を，
H, O には 6-31G** を用いた．
【結果・考察】１．静的水和モデル	  HF・B3LYP・MP2 レベルの計算を行ったが，いずれの結果も似通ったデータを与えたので，ここでは HF レベルの計算結果について言及する．構造最適化計算の結果，Mg2+, Ca2+, Sr2+, Ba2+ の第一水和圏に安定に存在し得る水分子数の数は，それぞれ，
6, 7, 9, 9 と見積もられた。これは，従来報告されてきた実験値・理論値（6.0, 6.0-7.9,
7.3-10.3, 7.8-9.3 [4]）とよく対応している．同様に Ra2+ の第一水和圏に幾つまで水分子が幾つまで水和できるかを検討したところ，9 個の水分子が水和した構造まで安定構造が見出された．しかし一方で，水和自由エネルギーは n = 6-9 において -267.6,-267.5, -266.8, -265.1 kcal mol-1 と求められ，水和数 n に対する依存性をあまり見せず，
静的モデルの計算結果のみからの Ra2+ の水和数決定は困難であった．
	
  NPA の結果より，軽アルカリ土類金属イオン（Mg2+, Ca2+, Sr2+）は，水和数（配位結合数）に比例した電荷減少 2.01.8 を示すが，重アルカリ土類金属イオン（Ba2+,
Ra2+ ）は，水和数の増加に伴う電荷減少量が 2.01.95 と少なかった．図 1 に
LMO-EDA の結果を示す．全てのアルカリ土類金属イオンの水和物において，金属と水分子の全相互作用のうち最も支配的な成分は静電相互作用（ES）項であり，2 番目に大きな相互作用成分は
電荷移動＋分極相互作用項（POL）であった．軽ア
ルカリ土類金属イオンから重アルカリ土類金属イオンになるにつれ，POL
項の寄与が小さくなり，
水和が配位結合的イオ
ン結合的に変化する様子
図 1 各アルカリ土類金属 [AE(H2O)1]2+ (AE=Mg, Ca, Sr, Ba, Ra) のがみてとれる．
-1
LMO-EDA 結果（kcal・mol ）２．動的水和モデル
	  図 2 に，水和 Ra2+ 系のRa-Hwater・Ra-Owater 動径分布関数を示す．Ra-Hwater・Ra-Owaterの第一ピークは，それぞれ2.85 Å，3.45 Å である．また，
水和数は統計平均の結果，約8 と算出された．図 2 から
Ra の第一水和圏の第一ピー
クは 2.55 ~ 3.55 Å の約 1.0
Å にわたってブロードであり，
Ra 周りの水分子はゆるく結合し，Ra まわりを動いている様子が伺え，静的水和モデ
ルからの議論と矛盾しない．
Ra2+ は，水分子と静電相互作用的な力で結合し，ある程度自由度をもった第一水分圏を形成することがわかった．

図2

水和 Ra2+ の動径分布関数 g(r)と水和数 CN(FMO3-HF/MCPdzp-CV(Ra);6-31G**(O,H))【参考文献】[1] P. Su, et al., J. Chem. Phys. (2009) 131, 014102., [2] Y. Komeiji et al., Chem. Phys. Lett., (2003) 372, 342.,[3] H. Anjima et al., J. Comput. Chem. (2007) 28, 2424., [4] R. Caminiti et al., Chem. Phys. Lett. (1977) 47, 275., A.
Tongraar et al., Chem. Phys. Lett. (2001) 346, 485., M. I. Bernal-Uruchurtu et al., J. Chem. Phys. (1995) 103, 1588. [5] A.
Tongraar et al., J. Phys. Chem. A (1997) 101, 6299., C. F. Schwenk et al., Pure Appl. Chem. (2004) 76, 37., S. Ramos et al.,J. Chem. Phys. (2003) 118, 5542., D. M. Pfund et al., J. Phys. Chem. (1994) 98, 13102., E. Spohr et al., J. Phys. Chem.
(1988) 92, 6754, T. S. Hofer et al., J. Phys. Chem. B (2006) 110, 20409, P. D'Angelo et al., Phys. Rev. B (1996) 54,12129, T.
S. Hofer et al., Chem. Phys. (2005) 312, 81