4Pp001

閉じ込められた量子系：
多電子量子ドットの電子構造(東大院理・マックスプランク宇宙物理)○佐甲徳栄・Geerd HF Diercksen【序】空間的に「閉じ込められた」量子系の研究は, 単にアカデミックな興味に留まらず, 量子ドットのモデルとして, フラーレンやナノチューブの中に拘束された原子・分子のモデル,
あるいは高圧力下に置かれた原子・分子のモデルとして, 近年大きな関心を集めている [1].
これまで我々は, 調和ポテンシャル中に拘束された原子 [2, 3] および分子 [4] の振電構造を調
べ, ポテンシャルの形状・閉じ込めの強さに応じて, 電子構造がどのように変化するのかを明らかにした. 本研究では, 調和ポテンシャル中に閉じ込められた複数の電子 (調和振動子量子ドット) [5] に着目し, そのエネルギー構造および電子密度分布を調べた.

【計算方法】本研究で用いたハミルトニアンを以下に示す：
n
n
n 
1 
Hˆ = ∑ − 12 ∇ i2 + ∑ 12 (ω x2 x i2 + ω y2 y i2 + ω z2 z i2 ) + ∑ 
.
r
r
|
−
|
i =1
i =1
i> j 

i
j


[

]

[

]

(1)

このモデルは, 三次元非等方調和ポテンシャルの中に閉じ込められた複数の電子のハミルトニアンを表す. 非等方閉じ込めポテンシャルに由来する様々な形状の電子分布を適切に表現するために, 次式で定義される非等方ガウス型基底関数系 [2-5] を用いて計算を行った.
r r
a
(2)
χ ξar (r ) = x a x y y z a z exp(−ξ x x 2 − ξ y y 2 − ξ z z 2 ) .
2 電子系の計算には完全 CI 法, 多電子系 (n = 3-6) の計算には多参照 CI 法を用いた.

【 Hartree-Fock 軌 道 】 spherical 型 ( ω x =ω y = ω z ), prolate 型 ( ω x = ω y > ω z ), oblate 型( ω x = ω y < ω z ), の三種類ポテンシャル形状について計算を行った. 図１に, ( ω x , ω y , ω y ) =(0.1, 0.1, 0.01) の場合の Hartree-Fock 軌道の軌道
エネルギーおよび電子分布を示す. この例は,異方性が大きい oblate 型量子ドットに対応する . 各軌道は, 二つの量子数 υ (= υ x + υ y ),

υ z , および z 軸に対する軸対称性のラベル ∆を指定して, [υ ,υ z ]∆ のように帰属される.
上下方向(z 軸方向)の強い閉じ込め, 水平方向の弱い閉じ込めを反映して, 図 1 に示すように低エネルギーの軌道は, 縦方向には節を持たず, z 軸を中心とする角度方向に節構造を持つことが分かる.
図 1. Oblate 型量子ドットの Hartree-Fock 軌道

【エネルギースペクトルと電子密度分布】図 1 で表される軌道に, 電子を 2 個および 4 個
つめた系, すなわち, 2 電子量子ドットおよび 4 電子量子ドットの, 一重項状態のエネルギー
スペクトルと電子密度分布を, それぞれ図 2, 3 に示す. これらの図が示すように, 一重項の最低エネルギー状態は, 2 電子系, 4 電子系ともに 11 Σ +g である。一方, 対応する電子密度分布は, 2 電子系は中心部に分布を持つ「円盤形状」をとるのに対し, 4 電子系では中央部の電子密度が極めて小さい「ドーナツ形状」を持つことが分かった. この形状の違いの由来を理解するために, CI 波動関数の主な配置を調べた. その結果, 2 電子系の場合には, 最低軌道に電子が 2 個入った ([0,0]σ g ) 2 が主配置であることが分かった. すなわち, 11 Σ +g 状態の電子密度分
布が示す円盤形状は, 図 1 に表される [0,0]σ g 軌道の形状を反映したものであると解釈される. 一方, 4 電子系の場合の主配置は, 最低軌道ではなく, 縮重した第一励起軌道に 4 個の電子が入った ([1,0]π u , x ) 2 ([1,0]π u , y ) 2 であることが分かった. 4 電子系 11 Σ +g 状態のドーナツ形状の電子分布は, 図 1 が表す中心部に節を持つ [1,0]π u 軌道の形状とコンシステントである.
これらの結果は, 量子ドットにおいては, 通常の原子の電子構造とは異なり, 低エネルギー状態においても多電子励起配置が重要となることを示している。

図 2. Oblate 型 2 電子量子ドットの電子構造

図 3. Oblate 型 4 電子量子ドットの電子構造[1] W. Jaskolski, Phys. Rep. 271 (1996) 1-66.
[2] T. Sako and G.H.F. Diercksen, J. Phys. B: At. Mol. Opt. Phys. 36 (2003) 1433-1457.
[3] T. Sako and G.H.F. Diercksen, J. Phys. B: At. Mol. Opt. Phys. 36 (2003) 1681-1702.
[4] T. Sako, I. Cernusak and G.H.F. Diercksen, to be submitted.
[5] T. Sako and G.H.F. Diercksen, J. Phys.: Condensed Matter (2003), in press.

4Pp002

LPN identity の分子系への応用と解析(筑波大化)[序]° 下堂 靖代 ,守橋健二

ESR スペクトル測定により観測される超微細結合定数 (hfcc) は、原子核上のスピン密度を直接反映する物理量である。近年、hfcc は、密度汎関数理論、Kohn-Sham(KS) 法により予測される場合が多くなりつつある。通常、核上のスピン密度は、局所的な演算子 δ 関数により求められるが、SCF の手続きで必要となる交換相関ポテンシャルは、核近傍で正しく振る舞わないことが報告されている [1]。また、従来、基底関数として用いられるガウス型関数 (GTO) による展開では、核近傍での性質が正しく近似されないことはよく知られている。したがって、DFT 計算で hfcc を予測する場合、このような核 cusp 上での不都合を避けることはできない。Liu-Parr-Nagy(LPN) [2]は、局所的問題の不都合を避け、全空間を積分して、核上の電荷密度を求める方法を提案している。
しかし、LPN identity は、スピン密度や分子系に応用された例はない。ここでは、LPN identity
を、原子や分子系の核上のスピン密度計算に応用した。

[ 計算方法 ]LPN identity では、核上の電子密度を計算する 4 つの等価な恒等式が導かれている。そのうち、最も容易に計算が実行できる式に基づくと、核上のスピン密度 ρs (0) は、以下のようになる。

ρs (0) = −ρs (r) =

1
4π

n
X

Z

r · ∇ρs (r)
dr
r3

(1)

Pijs χi (r)χj (r)ij

ここで、ρs および ∇ρs は、それぞれスピン密度関数と電子座標 r に関するその微分である。Pijs はスピン密度行列、χj は原子基底関数である。さらに、分子系における原子 N の核上のスピン密度

ρs (RN ) について、以下の式を導いた。分子系では、(1) 式は多中心積分となる。
1 X
ρ (RN ) = −
4π A
s

Z

wA (r)

X(r − RN ) · ∇ρs (r)dr=ρsA (RN )|r − RN |3
A(2)

ここでは、Becke の空間分割法を用い、一中心積分の和で表した。wA (r) は、原子中心 A のもつ
ファジーセル内での r における重み関数である。RN は原子 N の核座標である。

(1)、(2) 式を用いて、第一周期原子 (C, N, O, F) および開殻分子系 (CH, NH, OH) を対象とし、各核上のスピン密度計算を行った。スピン密度関数は、UHF 法および基本的な交換相関汎関
数 (Slater, SVWN, Becke, BLYP, B0-LYP) を用いた KS 法により求めた。基底関数は、原子系には、6-31G、Chipman の [5s2p] および [6† s4† p]、Dunning の [5s3p]TZ、Huzinaga の (11s7p) を、分子系については、Huzinaga の (11s7p/8s) を用いた。また、分子系の構造パラメータには、実験値を採用した。

[ 結果 ]

実際に数値計算を行ったところ、第 1 周期原子の ρs (0) および開殻分子系の ρs (RN ) は、それぞれ δ 関数演算子の期待値に等しい結果であった。なお、原子および分子系ともに、各核上の LPN 電荷密度もまた、δ 関数の期待値に等しい結果であった。したがって、LPN identity には、

δ 関数の値を改善するような効果は期待できない。しかし、LPN identity を、δ 関数演算子を作用させて求めた局所的な物理量を、非局所的な物理量として変換する方法とみなすことにより、以下のように従来の δ 関数演算子によるスピン密度計算では得られない情報を得ることができた。
ここでは、簡単のため、N 原子を例にとって説明をする。Table 1 には、(1) 式による ρs (0) 計算
結果を示した。Figure 1 は、N 原子の内殻領域について、核からの距離 r に対し (1) 式の被積分関数をプロファイルしたものである。原子のような球対称系では、(1) 式は以下のようにさらに簡単になる。

ρs (0) = −Z∂ρs (r)dr∂r(3)

LPN identity では、ρs (0) は全空間のスピン密度勾配と直接結びつく物理量であることがわかる。
局所スピン密度近似 (LSDA) による Slater、SVWN を除くと、核近傍でのプロファイルは正の領域から負の領域へと急激に変化している。UHF のプロファイルは、正の領域は最も高く、負の領域は最も浅くなり、(1) 式による ρs (0) は実験値を 2 倍近くも過大評価している。逆に、Slater、

SVWN のプロファイルは、核近傍で負の領域に留まっており、ρs (0) は負と、定性的にも誤った結果が得られる。LSDA では、この核近傍でのスピン密度の急激な変化を記述できていない。このことが、LSDA に定性的に誤った結果を与える原因となっている。交換項に GGA を加えること
により、この定性的な誤りは解消されるものの、BLYP でも実験値の半分程度の値しか得られな
い。Table 1 と Figure 1 を比較すると、正の領域でのピークが高く、負の領域での谷が浅いほど、

ρs (0) は大きくなっている。このことから、
ρs (0) の大きさは、主に、この核近傍でのスピン密度勾配の大きさに支配されていることがわかった。
さらに、このように、LPN identity を用いることにより、計算方法の違いによる ρs (0) および ρs (RN ) 計算値の相違を可視的に示すことが可能となり、また、これによりスピン密度値の相違の原因について明らかにすることができた。

Figure 1 : Integrand of Eq.(1) for the N atom

Table 1 : ρs (0) of the N atom obtained from calculating Eq.(1) with the (11s7p) basis set
pure-DFTsρ (0)hybrid-DFT

UHF

Slater

SVWN

SplusBBLYP

B0-LYP

exptl

0.1810-0.0550-0.01960.02120.06960.10480.09705

References[1] C. J. Umrigar, X. Gonze, Phys. Rev. A 50, 3827 (1994).
[2] S. Liu, R. G. Parr, A. Nagy, Phys. Rev. A 52, 2645 (1995).

4Pp003

デンドリマーの分極率の理論的研究(豊橋技術科学大学) ○岡田良平・関野秀男【はじめに】
ナノ分子をはじめとする巨大分子系の電荷ダイ
ナミズムのシミュレーションはますます重要となってきている。
図 1、図 2 のようなアゾデンドリマーの光学物性は非常に特異的であり、特に世代 4(L4)以上のアゾ
デンドリマーが示す光学物性の非線形性などに興味が持たれている。
しかし、そうした興味ある物性が発現する大きな世代のアゾデンドリマーは巨大分子であり、従来の図 1. trance-L3

量子科学計算で扱うのはかなりの計算労力を要する。また、アゾデンドリマーは柔軟な構造を持ち、その柔軟構造が特異的物性発現の鍵をにぎっている可能性があるが、すべての考えられる構造を最適化したり、その構造に対応する物性値をすべて理論算定するのは、現在のところ不可能である。
そこで、今回はあらかじめ最適化された代表的な構造のアゾベンゼン、世代 2(L2)、世代 3(L3)、世代 4 のアゾデンドリマーについて、その分極率の算
出をし、結果を分析した。

図 2. trance-L4【計算】
トランス型、シス型アゾデンドリマーの分極率の非経験的量子科学計算を行った。
Hartree-Fock 計算は Gaussian98[1]を用い、STO-3G と 6-31G**の基底関数によって行った。さらに、VWN 密度汎関数を用いた DFT 計算を、プログラムパッケージ ADF[2]を用いて、Slater 型基底関数(DZP)によって行った。
算出結果を表 1、表 2 に示す。(表の空欄は現在計算中)
表 1. トランス型アゾデンドリマーの分極率
trans-Azobenzene
trans-L2
trans-L3
trans-L4

ADF
Gaussian98
VWN/DZP
HF/6-31**
HF/STO-3G
180.01144.5592.62
690.16544.96348.551612.901308.24840.551749.16

表 2. シス型アゾデンドリマーの分極率
cis-Azobenzene
cis-L2
cis-L3cis-L4ADF
Gaussian98
VWN/DZP
HF/6-31**
HF/STO-3G
148.24125.0479.80
672.89546.80353.791559.901269.21821.151665.93【結果】
トランス型、シス型アゾデンドリマーの分極率の値をそれぞれ比較すると、いずれの世代もトランス型のほうが大きな値となっているが、世代が大きくなるとその絶対値が大きくなり、トランス-シス間の差の割合は小さくなることがわかる。
Gaussian98 による、HF/6-31G**、HF/STO-3G の算出値から、大きなサイズのア
ゾデンドリマーでも基底関数が分極率の算定にとって大きな影響を与えることがわかる。
密度汎関数はバンドキャップを過小評価し、巨大なポリエンなどの分極率計算で破綻することが報告されている[3]が、本研究では Hartree-Fock 計算とほぼ同様のサイズ依存性を示している。アゾデンドリマーの共役分子集合体としての特徴が原因であると考える。
前述したように、世代 4 以上のアゾデンドリマーが示す非線形性などに興味があり、
現在計算中の結果も含めて、その物性について解析している。
謝辞
本研究は、科学技術振興事業団の計算科学技術活用型特定研究開発推進事業(ACT-JST)(研究開発課題「DNA のナノ領域ダイナミクスの第一原理的解析」)の援助を受けて行われました。
参考文献[1] M.J Frisch et al., Gaussian98(Revision A.1) Gaussian Inc., Pittsburgh, PA,
1998[2] ADF2002.03, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam,
The Netherlands, http://www.scm.com[3] B.Champagne, E.A.Perpète, S.J.A.van Gisbergen, E.J.Baerends, J.G.Snijders,C.Soubra-Ghaoui, K.Robins and B.Kirtman, J.Chem.Phys, 109 (1998) 10489

4Pp004

強光子場中における分子のスピン状態がイオン化に及ぼす影響に関する理論的研究（東大院理，東大院工）

○鈴木

孝宗，河田

功，山下

晃一

【序】
高出力（10 14 W/cm 2 ）レーザーパルスは、それによる高精度物質制御の可能性をもたらすものである。高密度パルス光を原子に照射すると、Coulomb ポテンシャルが歪み、電子がトンネルイオン化する。イオン化した電子は光の電場から多光子分のエネルギーを得て再び核に戻る。この時、獲得したエネルギーを放出する形で、閾値イオン化（above-threshold ionization: ATI）や高次高調波発生（high-order harmonic generation:
HOHG）といった非線形現象が生じる。
分子の場合は、イオン化確率が平衡核間距離よりも長い核間距離Ｒ ｃ において最大となる enhanced ionization という現象が理論、実験双方の研究から見いだされている。このことは、核の運動に起因する新たな自由度が、分子におけるイオン化ダイナミクスにおいて重要であることを意味している。
しかし、分子におけるスピン状態がイオン化に及ぼす影響に関しては、これまでほとん
ど議論されてこなかった。そこで、本研究では核の運動のみならず、スピン状態の違いが強光子場中におけるイオン化にどのような影響を及ぼすか理論計算によって調べることにした。なお、簡便のため本研究は二電子系において行った。
【計算方法】
本研究は、二電子が分子軸方向（ｘ軸方向と定義）にのみ動き、核は等間隔に固定した一次元モデルを用いて行った。一次元モデルでは、時間依存 Schrödinger 方程式は原子座標系（ｅ＝ħ＝ｍ ｅ ＝１）においてi∂Ψ ( x1 , x2 , t ) = Hˆ Ψ ( x1 , x2 , t )∂tと表せる。ここで、Hamiltonian はＡ・ｐ相互作用を考慮して

 p12
  p22ˆ
H =
+ A ( t ) ⋅ p1  + + A ( t ) ⋅ p2  + V ( x1 , x2 )
 2
  2

とした［１］。なおＡ(ｔ)はベクトルポテンシャルでA ( t ) = − ∫ ε ( t ′ ) dt ′t0で与えられる。ここでε(ｔ)はレーザー電場である。本研究においては、レーザー電場は
分子軸方向に偏向した直線偏光を想定し、波長 800nm の矩形パルスを用いた。数値計算
はｘ １ ，ｘ ２ 双方向に 512 個のグリッド点を持つ 179.2a.u.のグリッド空間を用いて行った。
時間依存 Schrödinger 方程式は高速 Fourier 変換法により数値的に解いた［２］。核からの Coulomb ポテンシャルが有効である internal zone 中を波動関数が運動し、absorbingzone（この領域では Coulomb ポテンシャルはほぼ０）に達したら mask function を用いて波動関数を吸収した。このことによりグリッド端における波動関数の反射を防ぐことができる。時間発展が終了したら internal zone 中に残っている波動関数の絶対値からイオン化確率を計算した。
【結果】
Ｈ ２ および対称型直線状Ｈ ３ ＋ イオンのイ100singlettriplet示す。Ｈ ２ と対称型直線状Ｈ ３ ＋ ではイオン化の ス ピン 依存 性 が大 きく 異 なる 。Ｈ ２ の場合、一重項状態では、イオン結合性配置の形成により平衡核間距離よりも長い核間距 離 5.0a.u.の 時 に イ オ ン 化 確 率 が 最 大 とIonization probability

オン化確率の核間距離依存性を図１・２に
80
60
40
20
0

なるが、三重項状態では、核間距離が伸び2

るにつれイオン化確率が単調に減少した。
三重項状態では、二電子反発によりイオン

468Internuclear distance R(a.u.)

10

化が進行するので、核間距離増加に伴う基レーザー強度 1.5×10 14 W/cm 2 の時における底状態のエネルギー減少が外部障壁を超え

図１

るイオン化を抑制するためである。

Ｈ ２ のイオン化確率

一方、対称型直線状Ｈ ３ ＋ イオンの場合は、スピン状態によらずイオン化確率は似たよ
100の時にイオン化確率が最大となった。

80

このように同じ二電子系においても、核の数の違いによってスピン状態がイオン化に及ぼす影響が異なるのは非常に興味深い。
これらの結果から、水素分子と対称型直線状Ｈ ３ ＋ イオンでは enhanced ionization の

Ionization probability

う な 挙 動 を 示 し 、 最 外 核 間 距 離 が 9.5a.u.

60
40
singlettriplet2004

機構が異なることがわかり、核およびスピ

5

6 7 8 9 10 11 12 13 14 15Internuclear distance R (a.u.)ンの効果が重要であることが示唆された。
今後は、断熱場状態を用いて、イオン化直前におけるレベルダイナミクスを調べ、ス

図２

レーザー強度 1.0×10 14 W/cm 2 の時における対称型直線状Ｈ ３ ＋ イオンのイオン化確率

ピン状態の影響についてさらなる検討を行う予定である。
【参考文献】［１］R.London: The Quantum Theory of Light 2nd.ed. (clarendon,Oxford,1983) Chap.5［２］M.D.Feit et al , J.Comput.Phys. 47,412(1982)

4Pp005

単一有機分子の電気伝導に関する理論的研究（東大院工、産総研計算科学）丸山仁志、浅井美博、○山下晃一【序】近年、分子デバイスの開発に向けて、分子の電気伝導の研究が実験、理論の両方の立場から多くなされている。単一分子の電気伝導特性の詳細を明らかにすることは、分子レベルのエレ
クトロニクスを実現するために非常に重要である。実験的には、Reed らによって金電極に挟まれたベンゼンジチオレートの単一分子の電気伝導が測定されている。そこで、多くの理論研究がこの系を扱っており、本研究でも、まずはこの系を扱う。
【理論モデル】ベンゼンジチオレートのような共役分子ではコヒーレントな電気伝導が起こると考えられ、理論的には下のランダウアーの公式が考案され、これまでに大きな成果を収めている。

(111)面(111)面

2e 2
G=
T
h

電極 q
電極 p
Ｇがコンダクタンス、Ｔが透過確率、2e2/h がコン
図 1. モデル系
ダクタンスの量子単位である。本研究では、Fig.1のような半無限金電極の(111)面に挟まれたベンゼン
ジチオレートを考える。
ランダウアーの公式を発展させると、次のような電流Ｉが電圧Ｖの関数として計算される式にI (V ) =2e +∞dET ( E ){ fp( E − µ p ) − fq( E − µq )}
h ∫−∞

なる。
ここで、Ef はフェルミエネルギー、fp, fq は電極 p,q のフェルミ分布関数でμp,μq がその化学ポテンシャルである。また、透過係数ＴはエネルギーＥの関数であり、これをグリーン関数法によっ
て求める。
†
]
T ( E ) = Tr[ ΓpGall ΓqGall

Gall = ( E − H mol − ∑ p − ∑ q ) −1
Ｇall は系全体（分子と両側の電極）のグリーン関数であり、Hmol は分子のハミルトニアンである。Σp,Σq は自己エネルギーであり、半無限の電極 p,q の影響を擬ポテンシャルとして取り入れるものである。グリーン関数と自己エネルギーは tight-binding 法によって求め、本研究では拡張ヒュッケル近似を用いる。
【計算結果】図２にフェルミエネルギーを実験値の−5.3eV として計算した場合の電流−
電圧特性、図３にエネルギーに対しての透過確率のプロットと透過確率に寄与する分子軌道を示した。図３右図は透過確率をエネルギー

図２．ベンゼンジチオレートの電流−電圧特性に対してプロットしたもので、矢印はフェルミエネルギーを−5.3eV、電圧を 3Ｖとしたときのおおよその積分する範囲である。これを用いて、計算した電流−電圧特性が図２である。定性的な形は実験結果をおおよそ再現できている。図３の真ん中の図はベンゼンジチオレートの両側に
Au 原子を 4 個ずつ結合させたものの分子軌道の中で、軌道が分子全体に非局在化して透過確率に寄与していると考えられるものを抜き出してある。なお、分子軌道の横に書いてある数値はその固有エネルギーである。また、図３左図はベンゼンジチオレートの分子軌道の中で、軌道の位相関係から真ん中の Au と相互作用した場合の分子軌道の起源となると考えられるものを抜き出
している。

Ｅ[eV]

Au4C6H4S2Au4
の分子軌道

-1.97
S

S

S

S

S

S

S

S

-3.09

-4.37

-3.56

S

S

-4.38

Ef＝−5.3eVS-10.75

-4.80
S

S

S

S

S

S

S

S

-12.22

-7.05

-7.16

Transmission Probability
図３．透過確率のエネルギー依存性と透過確率に寄与する分子軌道

今後はベンゼンジチオレートで用いた方法を下のような非対称な分子に適用し、整流作用の可能性を探る。

NH

S

S

NO2

4Pp006
Energy Density Analysis (EDA)の GAMESS への組み込みといくつかの応用(早大理工) ○河村芳海、馬場健、中井浩巳
1. 緒言
近年の電子状態理論の発展により、量子化学計算は生体分子など巨大分子に対しても高精度に物理量を求められるようになった。しかし、通常これらの計算から得られるエネルギーは、単一のスカラー量であり、エネルギー移動過程など分子内部でのエネルギー変化をみることは困難である。当研究室で提案した Energy Density Analysis (EDA)は、量子化学計算で得られる全エネルギーを構成原子に分割する解析方法であり、これまでさまざまな系に適用することにより、種々の現象、反応メカニズムの解明に有用であることを示してきた[1-4]。当発表では、EDA を量子化学
計算パッケージ GAMESS.US に組み込み、さまざまな系への適用例を示すと共に、エネルギー分割法などと組み合わせた結果について述べる。
2. 理論
EDA では、次の手順で各エネルギー成分を構成原子に分割する。Kintic energy 項 Ts は、原子軌
道(AO)に基づいた分割を行い、atomic energy density TSA を得る。

TS = ∑∑ Pµν χν tˆ χ µ = ∑∑ Pµν Tνµ = ∑ (PT )µµ ,
µ

ν

µ

ν

µ

TSA = ∑ (PT )µµµ∈ A

Nucleus-electron attraction 項、Coulomb interaction 項、Exact exchange 項も同様に分割される。これにより Hartree-Fock エネルギーに対する分割が達成される。
この分割は Mulliken population analysis(MPA)の類推に基づいているが、同様に Löwdin population analysis (LPA)に基づいて直交化基底を用いた分割を行うこともできる。すなわち、{()(

TSA (L&o&wdin) = ∑ (P' T')µµ = ∑ S1 / 2 PS 1 / 2 S −1 / 2 TS −1 / 2
µ∈ A

µ∈ A

)}

µµ

Kohn-Sham 法においては、exchange-correlation 項を分点の重みωg、partition function pA を用いて次のように分割する。

= ∑ ω g p A (rg )FXC (rg )
grid

E

A
XC

O

g

C

3. 応用
3.1 表面吸着過程に対する EDA

H1

固体表面に対する分子吸着過程は触媒反応や電極反応の素過程として興味深い。本研究では Si(100)表面に対する CO 分子吸着過程を EDA

Si 1

Si 2

Si 2

H2

H2
Si 3

を用いて解析した。表面のモデルとしては、クラスターの端を水素原子でキャップする Capped bond model を用い、Si9H12 クラスター(Fig. 1)を用いた。基底関数は 6-31G(d,p)、DFT 法の汎関数は B3LYP である。

Si 1 H 1

H3
Si 4
H4

H4

Fig. 1. Si9H12 クラスター.

また、HF 法については MPA 及び LPA に基づく分Table 1. Si9H12 cluster への CO 分子吸着

割をそれぞれ計算した。吸着による各原子のエネルに伴うエネルギー密度変化 (kcal/mol).

ギー密度変化を EDA を用いて計算した結果を table

HF
Mulliken Löwdin
1.9
-25.2
-0.9
1.3
-0.4
-1.9
-2.5
0.6
0.5
0.1
0.1
0.3
0.2
0.6
0.3
0.1
-19.2
11.7
16.3
31.8
-2.2
-2.2

1 に示す。この系に対しては、HF 法について、MPAに基づく分割と LPA に基づく分割で定性的挙動が異なっている。Löwdin に基づく分割は、表面第 1層の Si dimer が大きく安定化し、吸着子が不安定化するという挙動を与えている。DFT 法における分割では HF 法における Löwdin と定性的に似た挙動を示している。
3.2 相互作用エネルギーの局所的な値への分割

Si1
Si2
Si3
Si4
H1
H2
H3
H4
C
O
Total

DFT
Mulliken
-32.4
0.8
-2.7
-0.6
0.2
0.2
-0.1
0.2
40.0
15.6
-10.3

水素結合など弱い結合により形成される分子系の相互作用は重要な研究分野であり、実験、理論両面から長年研究がなされている。量子化学計算において相互作用エネルギーは、孤立系と相互作用系のエネルギーの差分をとることで計算される。その相互作用の詳細な解析のために、相互作用エネルギーを物理的意味のあるコンポーネントに分割する手法が Morokuma により提案されて以来[5]、いくつかの分割手法が提案されてきた。Reduced Variational Space SCF method (RVS)
は変分空間を適当に制限して SCF を行うことにより、相互作用エネルギーを静電及び交換(ESX)、
分極(POL)、電荷移動(CT)相互作用コンポーネントに分割する手法である[6]。EDA と RVS を組み合わせることにより、相互H3or OH=2.050 A

作用エネルギーを物理的意味づけされた局所的な値に分割することが可能となり、分子間の相互作用の有用な解析手段のひとつとなることが期待される。Table 2 に、水二量体の相互

O1

H2

O4 H6
H5

Fig. 2. 水二量体.

作 用 に EDA/RVS を 適 用 し た 結 果 を 示 す 。 計 算 レ ベ ル は
HF/6-311++G(2d,2p)、単量 Table 2. 水二量体に対する EDA/RVS の結果 (kcal/mol). HF/6-311++G(2d,2p)体の構造は実験値に固定
ESX
POL
CT
Total Full variational
し、二量体間の距離を最適化した。相互作用エネルギ
ーのうち重要な静電及び交換相互作用は水素結合している水素と酸素に由来するものであることが

O1
H2
H3
O4
H5
H6
Total

1.94
-1.75
-1.37
-3.56
1.08
1.08
-2.60

-5.74
6.52
-0.90
-0.74
0.07
0.08
-0.70

-4.32
4.22
-0.40
-0.62
0.21
0.20
-0.72

わかる。
[1] H. Nakai, Chem. Phys. Lett., 363, 73-79 (2002).
[2] H. Nakai, K. Sodeyama, Chem. Phys. Lett., 365, 203-210 (2002).
[3] Y. Kawamura, H. Nakai, Chem. Phys. Lett., 368, 673-679 (2003).
[4] H. Nakai, K. Sodeyama, J. Mol. Structure (THEOCHEM), in press.
[5] K. Morokuma, J. Chem. Phys., 55, 1236 (1971).
[6] W. J. Stevens, W. H. Fink, Chem. Phys. Lett., 139, 15-22 (1987).

-8.12
8.98
-2.67
-4.92
1.35
1.36
-4.02

-8.51
9.52
-2.78
-5.26
1.48
1.48
-4.07

4Pp007

溶液内化学反応における分子論と連続体論の比較（京大院工）佐藤啓文・榊茂好【緒言】溶液内で起こる化学反応を調べるために、近年種々の電子状態理論が提案されている。各々の理論の有効性に対する議論は多々なされているが、それらを相互に比較・検討することは必ずしも十分になされているわけではない。本発表では、代表的な連続誘電体モデルであるPCM法とRISM-SCF法を取り上げて、それぞれの方法により求めた溶媒和された分子の電子波動関数の性質やエネルギー等について比較を行う。また、最近提案したエネルギー分配法を用いて溶媒和による分極の詳細について調べた結果について報告する。
【計算方法】具体的な比較は、溶媒和理論のベンチマーク系として最も広く用いられているSN2反応
Cl ＋CH3Cl →ClCH3＋Cl
について行った。気相中（孤立分子系）において、二つのCl-C間の距離の差
R=r(Cl-C)-r(C-Cl)
を固定し、他のすべての自由度について構造最適化してポテンシャルエネルギーを計算した。ついでこれらの構造を用いて、PCMおよびRISM-SCF計算を行い、水中における平均力ポテンシャルを計算した。PCM法における静電相互作用の評価には
GAMESSに実装されている標準値を用い、反発−分散力はFloris-Tomasiの表式でLennard-Jones型相互作用を採用した。RISM-SCF法の計算においても、これらと同じLennard-Jonesポテンシャルのパラメータを用いた。基底関数として6-311G*
を選び、すべての計算はHartree-Fock法を用いた。
25.0
20.0
gas phase

Energy / kcal mol-1【結果と考察】図１に気相におけるポテ
ンシャルエネルギーおよび二つの方法で求めた水中における平均力ポテンシャルを示す。既によく知られているように、気相で存在するイオン−双極子会合体（|R|=1~2Å）は水中ではほぼ消滅し、反応障壁は高くなる。PCM法およびRISMSCF法で求めた障壁の高さは極めてよい一致を示した。また、これらの値はこれまでに報告されている分子シミュレー
ションや液体の積分方程式理論の結果とも極めてよく一致した。一方でイオン−双極子会合体に相当する位置で、両者のエネルギープロファイルが異なっており、PCMに比較してRISM-SCF法では僅
かな安定化が見られる。

15.0

RISM-SCF
PCM

10.0
5.0
0.0
-5.0
-10.0
-10

-5
0
5
Reaction Coordinate / Å

図１：気相中および水中における
SN2反応のエネルギープロファイル

10PCM法、RISM-SCF法いずれの場合も溶液系の全エネルギー（Etotal）は溶質分子の電子状態に相当するエネルギー（Esol）と溶媒和自由エネルギー（∆µ）からなる。
Etotal=Esol+∆µ.
孤立状態における電子状態エネルギー（Eisolated）を導入すると、溶媒和に伴う溶質分
子の電子状態の変化に起因するエネルギー（Ereorg）は、Ereorg=Eisolated-Esol.

-30.0

5

Ereorg

-35.0

0

∆µ

−5

-45.0

−10

-50.0

−15

-55.0

−20

-60.0

−25

∆µ

-40.0

Ereorgとなる。図２に二つの方法による、 R にそった E reorg および ∆ µを示す。 ∆ µは絶対値に大きな違いがあるものの、全体では非常に似たプロファイルを示している。
すなわち反応原系および生成系に相当す
る|R|が大きな領域ではCl の溶媒和に伴う強い安定化が見られ、気相中での遷移状態に相当するR~0Åの領域で安定化は最小となる。絶対値の違いは系統的であり、反応原系からの差で相対値を比較す
るとほぼ等価となり、R~0Å附近でのRISM-SCF法の立ち上がりが比較的急激であることを除けば、両者に本質的な違いは無いと考えてもよい。一方、EreorgはR の全域に渡って寄与は小さいが、明確

−30

-65.0
-10

-5
0
5
Reaction Coordinate / Å10

図２：溶媒和自由エネルギーおよび電子状態変化
（実線：RISM-SCF法・点線PCM：法）が見られる。特に気相の会合体に相当する領域で、RISM-SCF法に比較してPCM法
では電子状態および溶媒和の変化が弱く緩やかであり、両者の平均力ポテンシャルのプロファイルの差異の主な原因となっている。
3.0

図３には、反応に伴うCl周辺の水の動径分布関数の変化を示した。実線はRISM-SCF法によって求めたもの、破線はPCM法で用いているcavityの半径に相当する距離から分布関数を階段関数として示したものである。Clに対する水素結合は、解離状態ではなく会合
体（R~1.5Å）で最も強くなっていることが分かり、これは先の議論とも一致する結果となっている。
当日はエネルギーの分割した結果についても議論する。

R=0.0Å

R=–10.0Å

2.5

Cl-O

2.0
1.5

Cl-H

1.0
0.5

cavitationelectrostatic dispersion-repulsion

0.0

R=10.0Å

R=1.5Å

2.5

2.0
1.5
1.0
0.5

0.0
0

1

2

3

4

5

6

7

0

1

2

3

4

5

6

7

8

図３：Cl周辺の動径分布関数
R<0ではClCH3に、R>0ではCl に相当している。

4Pp008
密度汎関数法による van der Waals 相互作用の解析（豊橋技術科学大学・工学部）○井上英計、栗田典之、関野秀男【はじめに】
密度汎関数(DFT)法は、分子の安定構造と電子状態を容易に且つ高精度に解析できるが、van der
Waals 相互作用のような分子間相互作用がどの程度記述できるかは、未だ確定した結論は得られていない[1-4]。本研究では、van der Waals 相互作用を定量的に記述できる密度汎関数を提案し、希ガス原子(He, Ne, Ar, Kr) 間の van der Waals 相互作用、van der Waals 相互作用の影響が大きいと考えられる DNA 塩基間の結合状態を、その密度汎関数を用いて解析し、その結果と伝統的非経験分子軌道法による結果を比較し、密度汎関数法で van der Waals 相互作用の解析がどの程度可能かを検討した。
【計算手法】
伝統的非経験分子軌道法 (MP2、MP3、MP4 及び CCSD 法)、及び DFT 法（PW91、Modified-PW91
等）により、希ガス原子間の相互作用エネルギーの距離依存性を調べた。更に交換相互作用の汎関数を調整し、van der Waals 相互作用を高精度に求めた。基底関数として、6-31G 基底に様々な分散基底を追加したものを用いた。また、DNA 塩基間のスタッキング相互作用の解析を行なった。
交換汎関数を調整した DFT 計算には NWChem プログラムパッケージ[5]を用いた。
【計算結果】
Table 1 に、各計算手法で得られた He-He の結合エネルギーと安定距離の値を示す。密度汎関数
法では、交換･相関ポテンシャル関数として、Perdew-Wang91 を用いた場合は、結合エネルギーが実験値の 10 倍にも過大評価されるが、PW91 を改良した Modified PW91 を用いるとこの過大評価は改善される。また、交換汎関数を調整すると、Modified PW91 よりも実験値に近い結合エネル
ギーを得ることができる。Figure 1 はシトシン２量体の結合エネルギーである。当日のポスターでは、Ne, Ar, Kr の結果と、DNA 塩基対の水素結合エネルギーとチミン２量体のスタッキング相互作用の結果ついても発表する予定である。
本研究は、科学技術振興事業団の計算科学技術活用型特定研究開発推進事業(ACT‑JST)(研究開
発課題「DNA のナノ領域ダイナミクスの第一原理的解析」)の援助を受けて行われた。

Energy (kcal/mol)Table 1 Binding energies and equilibium distances of He dimer
実験値 CCSD
Modified 調整交換
MP4 MP3 MP2 PW91
汎関数
[6]
(T)
PW91
Binding energy (kcal/mol) 0.021 0.021 0.018 0.017 0.012 0.210 0.056
0.054
Equilibrium distance (Å) 2.97
3.0
3.0
3.0
3.1
2.7
3.0
3.0
4.0
3.0
2.0
1.0
0.0
-1.0
-2.0
-3.0
-4.0
-5.0

xpw_d=3.4
xpw_d=4.0
xpw_d=4.4

2.8

3.0

3.2

3.4

3.6

3.8

4.0

4.2

4.4

4.6

Distance (Å)
4.8

Figure 1 Potential energy curve of Cytosine-Cytosine【参考文献】 [1] N. Kurita et al., Chem. Phys. Lett., 313 (1999) 693. [2] J. Sponer et al., J. Comput. Chem., 17 (1996) 841.
[3] D.C. Patton and M. R. Pederson, Int. J. Quant. Chem., 69 (1998) 619. [4] S. Tsuzuki and H. P. Luthi, J. Chem. Phys.,
114 (2001) 3949. [5] High Performance Computational Chemistry Group, NWChem, A Computational Chemistry Package
for Parallel Computers, Version 4.0.1 (2001), Pacific Northwest National Laboratory, Richland, Washington 99352, USA. [6]J. F. Olgivie and F. Y. H. Wang, J. Mol. Struct., 273 (1992) 277; 291 (1993) 313.

4Pp009
FMO-DFT による生体高分子に関する理論的研究（豊橋技科大 1、科技団 2、東芝 3）○杉木真一郎 1,2、栗田典之 1、関野秀男 1、田中茂典 3【序】
非周期的な巨大系を第一原理計算と同様の精度で計算することのできる Fragment MO（FMO）
法は、主として Hartree-Fock(HF)レベルでその有用性が検証されてきた。しかしながら生体物質
内で重要な役割を果たす水素結合などは、電子相関を考慮することが必要不可欠であり、FMO 法もそのような効果を考慮して計算する必要がある。今回、電子相関を考慮した電子状態計算法の一つである密度汎関数法を基にした FMO 法（FMO-DFT）を開発したので、これを用いて転写
活性タンパク質に関する研究を行う。
転写活性タンパク質の一つであるカタボライト活性化タンパク質 CAP（catabolite activator
protein）は、環状 AMP（cAMP）と結合し、CAP-cAMP 複合体となり、DNA の特定部位に結合し RNA ポリメラーゼが DNA に結合するのを助け、lac オペロンの転写を促進する。本研究では、CAP と cAMP の複合体の結合状態を研究することにより、CAP-cAMP 複合体がどのようにDNA と結合し、転写を促進するかを解明する足がかりとする。一方、cAMP を cGMP に変えた
場合、cGMP-CAP 複合体の構造が cAMP-CAP 複合体と比べて変化し、DNA との結合能が変化することが実験的に調べられている。そこで、cGMP-CAP 複合体についても、cAMP-CAP 複合体と同様に計算を行い、両者の比較をすることによって、その特性について理論的側面から検証
を行う。
【理論】
FMO は全系の計算を行わず部分空間に切り分け、それぞれ独立に計算する方法である。こうする
ことによって少ない計算労力で巨大な系を計算することが可能になる。
簡単に FMO の計算方法について HF の場合を例に用いて概略を示す。最初に分子を N 個のフラグメントに分割する。次にフラグメントをそれぞれ独立に計算する。フラグメントを計算する
際に、計算するフラグメント以外の（N-1）個のフラグメントからの静電ポテンシャルをフラグ
メントに足しこむ。

fI =−

all
N
ZA
1 2 nI / 2
ρ (r ' )∇ + ∑ {2 J a − K a } − ∑+ ∑ ∫ dr '+ ∑ Bc θ c θ c2| r1 − r ' |A | r1 − r A |K (≠ I )ca(1)

ここで、nI はフラグメント I の電子数、K は I 以外のフラグメント、Bc は十分大きい正の値、θ
cは局在化された軌道であることを表す。周囲からの静電ポテンシャルによってフラグメントの電子密度行列は変化するのでこれが収束するまで繰り返す。
次にフラグメントペアの計算を行う。フラグメントペアの計算の際には、フラグメントの計算の場合と同様に周りの（N-2）個のフラグメントからの静電ポテンシャルを足しこむ。

n IJ / 2
allNZAρ (r ' )1
f IJ = − ∇ 2 + ∑ {2 J a − K a } − ∑+ ∑ ∫ dr '+ ∑ Bc θ c θ c| r1 − r ' | c
2
A | r1 − r A |K (≠I ,J )a(2)

ここで、nIJ はフラグメントペア IJ の電子数である。フラグメントペアの計算は一回のみ行う。
式(1)、式(2)を用いてそれぞれフラグメント、フラグメントペアのエネルギーEI,EIJ を求める。その後に系のトータルエネルギーE を下式を用いて求める。
N

N

J >I

I

E = ∑ E IJ − ( N − 2)∑ E I(3)

次に FMO 法の計算結果より得られる有効ペア相互作用について、以下に簡単に示す。
Fragment および Fragment-Pair のエネルギーを Fragment 自身のエネルギーと周りからの静電ポテンシャルに分けると次のように書ける。
N

E I = E I' + ∑ V I ← K

(4)

K ≠I

N

E IJ = E IJ' + ∑ VIJ ← K(5)K ≠I ,J

ここで、E’I,E’IJ はそれぞれ Fragment および Fragment-Pair 自体のエネルギーであり、V は周りからの静電ポテンシャルになる。式(4),(5)を全エネルギーの表式(3)に代入し、整理すると、
N

N

I

IJ

E total = ∑ E I' + ∑ ∆E IJ(6)となる。∆EIJ は有効ペア相互作用と呼ばれるものであり、
N∆E IJ = ∆E IJ' + ∑ ∆VIJ' ← K(7)∆E IJ' = E IJ' − E I' − E J'

(8)

∆VIJ' ← K = VIJ ← K − VI ← K − V J ← K(9)K ≠I ,Jとなる。∆EIJ の値を各フラグメント間で比較することにより、フラグメント間の相互作用を定性的に見積もることができる。今回の研究において、CAP と cAMP（cGMP）の結合状態を見積もる方法として、FMO 法の結果より得られる有効ペア相互作用を用いる。有効ペア相互作用を見る
ことによって、cAMP と CAP の残基との相互作用を定性的に見ることができる。
【方法】
Protein Data Bank（PDB）に登録されている CAP-cAMP の構造（PDB コード：1G6N）を入手し、cAMP とその周辺のタンパク質を切り出す。そして、タンパク質側は２残基を１フラグメ
ントとし、cAMP および cGMP はそれ自体を 1 フラグメントとして FMO 計算を行う。そして、
その FMO の結果から得られた有効ペア相互作用の結果を用いて、cAMP 及び cGMP と CAP との特異的な相互作用について論じる。詳細については当日報告する。また結晶水の影響などについても議論を行う予定である。
【謝辞】
本研究は、科学技術振興事業団・計算科学技術活用型特定研究開発推進事業の研究課題「DNA の
ナノ領域ダイナミクスの第一原理的解析」の援助を受けて行われた。

4Pp010
局在化軌道を用いた電子相関理論の開発と応用（京大院工1、東大院工2） ○中尾 嘉秀1、榊 茂好1、平尾 公彦2【序】 局在化軌道を用いた電子状態理論は動的電子相関を効率的に取り込むことができるため、大規模な分子系に適用できる理論として興味が持たれている。これは動的電子相関が本来、電子同士の近接という局所的な現象から生じるため、分子全体に拡がる軌道よりも局在化軌道の方が記述に適しているためである。実際の計算コストを考えると、canonical 軌道を用いた場合では AO
積分から MO 積分への積分変換に基底関数サイズの 5 乗に依存して計算コストが増大する。局在化軌道を用いた場合は低いコストで積分変換が可能となるが、エネルギーを見積もるためには大きな計算コストを必要とする。そこで、本研究では大規模な分子系への電子状態理論の適用を目指して、局在化軌道を用いた MP2 法とそのプログラムの開発を行うことを目的とする。
【理論】 occupied 軌道は Boys localization 法により得られる orthogonal 軌道を、virtual 軌道として(1)式で示されるように、AO を virtual 空間に射影して得られる non-orthogonal 軌道を用いる。
φ= (1 − P ) χ AO

(1)

積分変換は Fig. 1 にあるように、占有軌道 i に近接する仮想軌道の組{a}に対してのみ、また占有
軌道 i と j の組に対しては全て変換を行う。

(S 2 )ia ,i′a′ = δ ii′ S aa′（2） (F2 )ia ,i′a′ = Fii′ S aa′ − δ ii′ Faa′（3）(2)式で定義される S2 行列を対角化するユニタリー行列 US を求めて、固有値が小さいベクトルを除去し、規格化する行列 XS を求める。次に、(3)式で定義される F2' (=XS†F2XS)を対角化する行列
XF を求めて、固有値が小さいベクトルを除去する。以上の操作により、仮想軌道が持つ線型従属の問題が除去できる。X(=XSXF), f2 (=X†F2X), V' (=X†VX)を用いて MP2 の 2 次のエネルギーE2 は(4)式のようになる。
′ Viajb′
E 2 = ∑ Viajb{(f ) + (f ) }
2 ia

2

jb

(4)

iajb

【計算と結果】 Table １に cis-retinal, trans-retinal, hexahelicene, α-carotene, β-carotene, taxol の局
在化軌道を用いた MP2 法の 2 次の摂動エネルギーを記した。基底関数は 6-31G を用い、virtual 軌
道は occupied 軌道から 2.0Å 離れた軌道までを含めた。その結果、分子のサイズに依存せず、約 25 個のvirtual 軌道 a

virtual 軌道 b

occupied 軌道i

virtual 軌道を用いるだけで電子相関の 約 98%を見積もる。計算時間はtaxol で canonical 軌道を用いた場合の 10 分の１以下となり、非常に効率

核
occupied 軌道ｊFig. 1.

Occupied orbitals and virtual orbitals.

よく電子相関を取り込むことができる。また、本研究で開発されたプログラムは並列化されており、Table 2 に Taxol を 6-31G(660 基底)、6-311G**(1484 基底)を基底関数として用いた場合の計算時間を示した。積分計算と初めの 2 インデックスの積分変換は完全に並列化されており、使用した CPU 分の性能を発揮している。残りの 2 インデックスの積分変換とエネルギーの評価の際には、ノード間のデータ転送がボトルネックとなって効率が低下する。8 CPU を使った計算では約 6 分の 1 の時間で計算可能である。
Table I.

The correlation energies by the localized MP2 method (basis 6-31G).

Molecule

Occ.

K

E(HF)

E2(canonical)

Dom.

E2(localized)

cis-retinal

78

245

-848.11443-1.92143

23.2-1.88605(98.2%)

trans-retinal

78

245

-848.12502-1.91919

23.2-1.88501(98.2%)

hexahelicene

86

266

-993.58719-2.2522726.0-2.18475(97.0%)

α-carotene

148

472

-1546.60474-3.59358

23.8-3.52361(98.1%)

β-carotene

148

472

-1546.61685-3.59522

23.8-3.52260(98.0%)

taxol

226

660

-2910.69150-6.00386

24.2-5.87186(97.8%)

Occ.: Number of occupied orbitals; K: Number of basis functions; Dom.: Number of domain orbitals per
occupied orbital; E2: The second order correlation energy of MP2.

Table II.

Timing in seconds of the parallelized local MP2 calculation (wall clock time). The values in theparentheses denote the number of effective CPUs.
CPUs

Int+Trf12a

Trf34bCalc Ec

Total

Taxol (6-31G basis set, 660 basis functions, 164 correlated electrons)
1

4,729.4

2,564.9

9,881.1

17,219.7

22,389.8 (2.0)

1,684.0 (1.5)

5,206.5 (1.9)

9,321.7 (1.8)

41,232.1 (3.8)

1,246.8 (2.1)

2,854.5 (3.5)

5,379.5 (3.2)

8613.3 (7.7)

1,087.7 (2.4)

1,862.7 (5.3)

3,626.8 (4.7)

Taxol (6-311G** basis set, 1,484 basis functions, 164 correlated electrons)
1

66,777.1

50,096.4

103,702.0

221,249.8

233,731.7 (2.0)

27,374.9 (1.8)

56,709.7 (1.8)

118,293.9 (1.9)

416,838.8 (4.0)

16,798.0 (3.0)

31,247.7 (3.3)

65,416.7 (3.4)

88,469.7 (7.9)

11,761.2 (4.3)

17,475.0 (5.9)

38,342.3 (5.8)

aCalculations of two-electron repulsion integrals and first half transformation.

b

Calculation of second half transformation.

c

Calculations of solving the T vector and evaluating the correlation energy.

r4Pp011直線型及び分岐型分子集合体におけるエキシトン移動（阪大院理 1・阪大院基礎工 2）○高畑昌弘 1・中野雅由 2・山田悟 1・岸亮平 2・新田友茂 2・山口兆 1【序】近年、励起エネルギ− 移動において、高分子や分子集合体の特異な構造との関連において興味深い現象が見られている。例えば、ケ− リ− ツリ− と呼ばれる構造を持ったフェニルアセンチレンデンドリマ− は、その分子の外縁部で吸収した光エネル
ギ− を中心へと集める光エネルギ− 収穫機能を持つことで知られている。これらの分子は、外縁部から中心へ向けてπ共役の鎖長が段階的に長くなり、各直鎖部が接合する部位で共役が切断されている特徴を持つことから、励起エネルギ− が高い状態では外縁部にエキシトン分布を多く持ち、低い状態では中心部に分布を多く持つというエ
キシトン状態を持つ。また、励起状態間における遷移は、エキシトン分布の重なりの大きな状態間で起こりやすいことから、外縁部に分布の多い状態から、重なりの大きな中間部分に分布の多い状態への遷移が起こり、最後中心に分布の多い状態への遷移が起こることが分かった。このような段階的な遷移の結果として、空間的には外縁部
から中間部分を経て中心へと励起エネルギ− が移動することが分かった。これまでの研究から、（ａ）部分的には非局在的な状態を形成しながらも、それらがお互いによく分離されていることと（ｂ）それらの非局在化した部位が、段階的なエネルギ− 移動を起こしうるように、すなわち、各分布が順次に重なりを持つように、配置されていることの二つの条件を満たすならば、ケ− リ− ツリ− 構造に限らず、方向性のある
エネルギ− 移動を起こしうる可能性があることを予測した。本研究では、その一つの例として、直線型分子集合体を取り上げる。
【方法】図１に示した二つの型の直線型分子集合体を考える。ここで、矢印はモノ
55 chromophores
1

2

3

10

(A)

……
10 a.u.10 a.u.
2
1

10

3

……

(B)
20 a.u. 20 a.u.

図１

二つの直線型分子集合体の構造。（Ａ）モノマ− 間の距離が一定（10 a.u.）の直線型分子集合体。（Ｂ）各セグメント間の距離（20 a.u.）をモノマ− 間の距離（10 a.u.）の２倍にした直線型分子集合体。

(A)-2, 3, 4, 5
2

4

2.0

3

5

1.0

2.0

2

4

3

5

1.0

0

2
Time [ps](A)-7, 8, 9, 10
7

4

0.0

9

3.0

0

2
Time [ps](B)-7, 8, 9, 10

8

10

-4

]

4.0

-4

]

0.0

(B)-2, 3, 4, 5

]

3.0
-4

-4

]

3.0

7

9

8

10

4

2.0
2.0
1.0
0.0

図２

0

2
Time [ps]

4

0.0
0

2
Time [ps]4直線型分子集合体モデル（Ａ）と（Ｂ）のエキシトンダイナミクスの結果。各セグメントに属するモノマ− の寄与を足し合わせたエキシトンのポピュレ− ションの変化。

マ− の遷移双極子モ− メントを表す。各モノマ− は２状態からなると仮定し、分子間には双極子− 双極子相互作用を考慮する。モデル（Ａ）はモノマ− 間の距離が一定（10a.u.）であり、各セグメントに分けて示したのは、モデル（Ｂ）との対照のためである。モデル（Ｂ）は各セグメント内のモノマ− 間距離は 10 a.u. であり、各セグメン
ト間の距離は 20 a.u. である。これら二つの分子集合体モデルに対して、電場印加後のエキシトンダイナミクスを、Born-Markov 近似の下で、電子− 格子相互作用を考慮したマスタ− 方程式を数値的に解くことによって実行する。
【結果】図２にエキシトンダイナミクスの結果を示す。ここではモデル（Ｂ）との対照のために、各セグメントに属するモノマ− の寄与を足し合わせたエキシトンのポピ
ュレ− ション変化を示す。この結果から、モデル（Ａ）において十分時間が経過した
後は、エキシトンの移動が見られないことが分かる。これに対して、モデル（Ｂ）においては、各セグメント間のエキシトンの移動が見られ、十分時間が経過した後は、
モノマ− が９個又は１０個属するセグメントのポピュレ− ションが大きくなっていることが分かる。これは、非局在化とその分離の効果から、非局在化のために励起エネルギ− が低いセグメントへの段階的なエキシトン移動が起こっていると考えられる。このモデルに対する相対緩和因子による詳細な解析及びこの発展的なモデルとして分岐型分子集合体モデルについて当日発表する。
【参考文献】[1] M. Takahata, M. Nakano, H. Fujita, K. Yamaguchi, Chem. Phys. Lett. 363 (2002) 422.
[2] M. Takahata, M. Nakano, K. Yamaguchi, J. Theor. Comp. Chem. Submitted.

4Pp012

尿素・アセトンへの陽電子吸着-配置間相互作用法を考慮した多成分分子軌道法によるアプローチ（横浜市大院総合理学）○立川 仁典【序】近年、陽電子を利用した陽電子分光法が、固体表面や格子欠陥を探るプローブとして注目を集めている。医療への応用も盛んとなり、陽電子放出断層法(Positron Emission Tomography,
PET)が癌診断に利用され、実用化レベルとなったのは記憶に新しい。このような技術革新に伴い、陽電子が中性分子に捕獲され、安定な共鳴状態生成を示唆する例も報告されている[1]が、陽電子照射による電子・陽電子状態の変化や、構造緩和に関する知見は殆ど解明されていないのが現状である。一方、我々は今までに、電子だけでなく、プロトン、ミューオン、陽電子といった軽い粒子の量子効果を直接取り込むために、分子軌道の概念をそれら粒子にも拡張させた「多成分分子軌道法」を開発・展開してきた[2]。本手法は、Born-Oppenheimer 近似が成立しないような分子系や、陽電子化合物を解析するのに強力な手段となる。独立粒子近似から出発すると、電子相
関だけでなく、陽電子−電子相関といった異種粒子間の多体効果を考慮する必要がある。本発表
では、このような多体効果を取り込むための「多種粒子状態理論の CI 法への拡張[3]」、および本手法の「尿素･アセトン分子への陽電子吸着[4]」といった内容を報告する。
【理論】全ハミルトニアン(Ｈtot)を、電子部分(Ｈe)、陽電子部分(Ｈp)、電子−陽電子部分(Ｖep)に分割する。全波動関数(Ψtot)は、電子・陽電子各波動関数(Φel・Φpm)を用いて、Ψtot = ΣΣΦel Φpm
C(l,m)と展開する。CI 係数 C(l,m) を求めるには、電子部分の行列要素( ＜Φel| Ｈe |Φel’＞ )だけ
ではなく、陽電子( ＜Φpm|Ｈp|Φpm’＞ )、電子−陽電子間( ＜Φel Φpm|Ｖep|Φel’ Φpm’＞ )の行
列要素も計算する必要がある。本計算では、HF 法および CI 法により、尿素分子・アセトン分子に陽電子を付着させた系を考察した。電子には aug-cc-pvdz、陽電子には充分な diffuse 関数を考慮した[9s5p1d]基底関数を用いた。CI 計算の励起配置として、1 電子、1 陽電子、および 1 電子
一 1 陽電子、各励起配置を含めた。

(a)【結果】図１に HF/aug-cc-pvdz レベルで最適化した尿HHN素・アセトン各分子の構造を示す。それぞれ C2, C2v 対C

メント、エネルギーを示す。アミノ基とメチル基の違いH

子に陽電子を付着した計算を実行した。単体のエネルギ

H

(b)

H
H
H

C

ーから陽電子化合物のエネルギーを差し引くことにより、陽電子親和力(PA)の値を算出することができる。HF
および CI 法により求めた PA を表２に示す。全ての PAが正の値となっているが、これは陽電子が付着された状

Table 1. Energy and dipole moment of urea
and acetone molecules with HF/aug-cc-pvdz
Species
Dipole[D]Energy [hartree]
Urea
3.990
-224.0159205
Acetone3.263-191.9728030

O
-0.68

N -0.18により、尿素の双極子モーメントがアセトンよりも大きくなっているのがわかる。次に図１の構造の下で、各分

1.200

1.23

称性を有する。表１に最適化された構造での双極子モー

C

1.193

0.76

H
H

O

-0.60

C -0.24
H

Figure 1: Optimized structure of(a) urea and (b) acetone moleculeswith HF/aug-cc-pvdz.

Table 2. Energy and positron affinity (PA) for positronic urea and acetone
molecules with HF and CI
Species
Positron basis Wave function Energy [hartree] PA [kcal/mol]+[9s5p1d]
HF-224.0161437
0.14[Urea ; e ][11]-224.0164133
0.31[9s5p1d]CI[Acetone ; e+][9s5p1d][9s5p1d]

HF
CI[11]-191.9728377-191.9729683

0.02
0.10

態が安定に存在することを意味する。 CI 計算により、尿素で 0.17kcal/mol、アセトン で
0.08kcal/mol 程 PA が改良されているが、これは陽電子−電子相関が重要な寄与を示すためである。図２にそれぞれの分子の電子分布、陽電子分布を示す。電子分布がカルボキシル基の酸素側に偏っているため、陽電子分布も酸素側に大きく偏っているのがわかる。つまり陽電子の付着機
構は、ポジロトニウム（電子と陽電子の複合体）生成というよりも、むしろ分子から生じる双極子場によって陽電子が捕獲される寄与の方が大きいことがわかった。実際、双極子モーメントの大きな尿素分子の方が、陽電子分布が局在化しており、陽電子親和力（陽電子状態）と双極子モーメントの間に相関を見出すことができる。本手法により得られる結果は、陽電子を含めた様々な実験や医療レベルでの先端治療に対して、多くの基盤的情報を与え得るものと期待される。
16

16

8

8

y (au)y (au)

(a)

0

-8

-8

-16
-4

0

4

12

x (au)

20

28

-16
-4

16

16

8

8

4

12

20

28

12

20

28

x (au)y (au)y (au)

(b)

0

-8

-16
-4

0

-8

4

12

x (au)

20

28

-16
-4

4

x (au)

Figure 2: Electronic and positronic density of (a) urea and (b) acetone molecules.

[1] Sueoka and M. Kimura, Phys. Rev. A, accepted (2001). M. Tachikawa, I. Shimamura, R. J.
Buenker, and M. Kimura, New directions in Antimatter Chemistry and Physics, Chapter 23,
437, Kluwer Academic Publishers. 2001. [2] M. Tachikawa, K. Mori, K. Suzuki, and K. Iguchi,
Int. J. Quantum Chem. 70, 491 (1998). [3] M. Tachikawa, Chem. Phys. Lett. 360, 494-500(2002). [4] M. Tacikawa, R. J. Buenker, and M. Kimura, J. Chem. Phys. 119, xxxx (2003).

4Pp013
フラグメント分子軌道による高次の物性量算定(豊橋技術科学大学大学院

知識情報工学系)○澤田寛明、杉木真一郎、関野秀男【はじめに】
現在、生体高分子やナノスケールの巨大分子を対象とした電子状態計算が盛んに行われている。このような巨大系を対象とした場合、近似なしの Ab-initio 計算は困難である。そして、大規模系に適用できる精度の高い Ab-initio 計算の近似手法として、
フラグメント分子軌道(FMO)法[1]が開発された。
ＦＭＯ法では、全エネルギーはフラグメントのエネルギーとフラグメント対のエネ
ルギーの重ねあわせで計算され、その全エネルギーの計算誤差は近似なしの Ab-initio
計算と比較して数 kcal/mol 程度であることが知られている。
ＦＭＯ法を用いて水、アンモニアの 3 量体についての１次、2 次の物性量を算定したところ、近似なしの Ab-initio 計算と比較し、1 次では数%、2 次では 0.3%程度の誤差で物性量が得られた[2]。これにより、ＦＭＯ法は占有状態のみならず、非占有状態についても高い精度で算出できることがわかった。
本研究では、さらに高次の物性量についてＦＭＯの精度を検討する。
【理論】
ＦＭＯ法における密度行列は
ρ FMO = ∑ ρ I ⊕ ∑ ∆ρ IJ
I

I <J∆ρ IJ = ρ IJ − ρ I − ρ Jで求める。ここで ρ I , ρ J , ρ IJ はフラグメント I、J とフラグメントペア(I、J)の密度行列である。
高次の分極率は系の外場に対する高次の応答として定義される。

µ = µ 0 + αε +11
βε 2 + γε 3 + L
2
6

ε : 外部電場

ここでは、上のようにして得られた密度行列の外場に関する微分から高次の分極率
β , γ を算出する。
【方法】
対象分子は図 1 に示すような、水(W)とアンモニア(A)の 3 量体および、図 2 に示すようなシトシン(C)とグアニン(G)の 1 本鎖 DNA を対象とした。計算は Ab-initio 計算
プログラム NWChem[3]にＦＭＯを実装したものを利用し、基底関数は STO-3G を用いた。
【結果および考察】
表 1 に水、アンモニアの 3 量体について、表 2 にシトシン、グアニンの 1 本鎖 DNA

の双極子モーメント( µ )、分極率( α )および超分極率( β , γ )の結果を示す。
表 1 より、W-A-W のような分子集合体では 2 次の超分極率まで数%以下という誤差の範囲で算出できていることがわかる。
3 次の超分極率でも主成分の有効桁数(2 桁）を考えると、有効な近似であることがわかった。

図1：水（W）-アンモニア（A）-水（W）の構造

表 1：W-A-W 構造の結果
HF-FMOHF

図2：シトシン（C）-シトシン（C）-グアニン（G）の1本鎖DNAの構造
表 2：C-C-G の 1 本鎖 DNA 構造の結果

error

HF-FMO

HF

error

µx

-0.50696 -0.50948-0.4942%µx-11.833-11.9413-0.9071%µy-0.26265 -0.26422-0.5936%µy0.489398

0.480895

1.7681%µz
α xx

0.06922 0.068808

0.5989%2.406886

2.452537

1.8614%10.03678 10.00166

0.3511%µz
α xx

325.0231

324.8534

0.0522%

α yy

9.003374 9.005107

0.0192%

α yy234.0709

233.5270.2329%

α zz
β xxx

9.51107 9.502679

0.0883%179.6191

179.6536

0.0192%5.757286 5.383839

6.9364%

α zz
β xxx

56.54601-68.0302-183.1190%

β yyy-11.2028 -10.8676-3.0850%

β yyy-225.093-100.42-124.1518%

β zzz
γ xxxx-4.50171 -4.48705-0.3268%-60.0485-35.9159-67.1917%

340.1226 300.4175

13.2167%

β zzz
γ xxxx25502.2320252.63

25.9206%

γ yyyy128.6367 103.5687

24.2042%γ yyyy-6018.58

7863.491

176.5383%

しかし、1 本鎖 DNA 構造などでは分極率の近似はかなり良いが、超分極率はうまく近似できていないことがわかった。これは、主鎖を通した電子移動のメカニズムが超分極率の算定に重要であるためと思われる。
【参考文献】
1：K.Kitaura,T.Sawai,T.Asada,T.Nakano,M.Uebayashi Chem.Phis.Lett.312(1999)319
2：(豊橋技科大) 澤田寛明、仙石康雄 、杉木真一郎、関野秀男
3 ： Richland ， Wachington ， Pacific
Version4.1(2002)

Northwest

理論化学討論会(2003)

National

Laboratory ， USA ．

4Pp014
QM/MM レプリカ交換を用いた最安定構造の探索(1 分子研･2JST/CREST･3 総研大) ○1, 2 川島雪生･1, 2, 3 岡本祐幸
【序】大きな分子クラスターには、局所安定構造が非常に多く存在し、通常の geometry
optimization algorithm では最安定構造を求めることが困難である。それに対して、レプリカ交換法 1-4 はポリペプチドやクラスターのような局所安定構造が多数存在する系の最安定構造を求めるのに非常に有用であることが示されている。
本研究では、このレプリカ交換法を分子内包クラスターに適用する。近年、有機合成などの分野において分子内包クラスターの研究は非常に盛んであり、X 線回折による構造解析が行われているが、水素原子については調べることが出来ない。そのため、内包されるもの（ゲスト）が水クラスターなどの場合、酸素分子の位置しかわからず、まわりのクラスター（ホスト）との相互作用などの議論が出来ない。
そこで、本研究は、クラスターの最安定構造の探索し、ホストとゲストの相互作用について調べることを目的とし、レプリカ交換法を水クラスターが内包される分子内包クラスターに適用する。

【計算手法】系が大きいため、すべての分子を ab-initio 分子軌道計算で扱うのは難しい。そこで、ゲストの水分子を Quantum Mechanics、ホストを Molecular Mechanics で扱う QM/MM法を用いる。ホストとゲストは 3Å以上はなれており、ホストを MM で扱っても差し支えない。ホストはゲストと比較し、非常に重いのでゲストの配置のみをモンテカルロ法にてサン
プリングする。徐冷法とレプリカ交換法を用いてシミュレーションを行う。本研究の計算には、UTChem5 を用いており、量子化学計算では INTEGRA6 module、シミュレーションでは、
S&D7 module をそれぞれ用いている。
【結果】当日は古典力学を用いた結果と QM/MM 法を用いた結果とを比較する。また、徐冷法による結果とレプリカ交換法による結果についても議論する。
[1] Hukushima K. and Nemoto K.: J. Phys. Soc. Jpn. 65 (1996), 1604
[2] Geyer C.J.: Keramidas E.M. (Ed.), Computing Science and Statistics: Proceedings of the 23rd
Symposium on the Interface, Interface Foundation, Fairfax Station (1991) 156[3] Swendsen R.H. and Wang J.-S.: Phys. Rev. Lett. 57 (1986) 2607
[4] Ishikawa Y., Sugita Y., Nishikawa T., Okamoto Y.: Chem. Phys. Lett. 333 (2001) 199.
[5] UTChem: Yanai, T.; Kamiya, M.; Kawashima, Y.; Nakajima, T.; Nakano, H.; Nakao. Y; Sekino, H.;
Paulovic, J.; Tsuneda, T.; Yanagisawa, S.; Hirao, K.
[6] INTEGRA: Emoto, M.; Iikura, H.; Inoue, H.; Kamiya, M.; Kawashima, Y.; Nakajima, T.; Nakayama,K.; Sekino, H.; Sorakubo, K.; Tawada, T.; Tsuneda, T.; Yagi, K.; Yanagisawa, S.; Yanai, T.
[7] S&D: Kawashima Y.

4Pp015

酵素モデル化合物の溶液中における物性の理論的研究(京大院理) ○東 雅大・加藤 重樹
遷移金属を含む酵素は、自然界で最も重要な役割を担う物質の一つである。これらの物
質では、活性部位に遷移金属を複数含むものが多く、その遷移金属間の相互作用が反応や物性にどのような影響を与えるか興味が持たれている。この分野の研究は古くから行われており様々なデータがあるものの、電子状態を理論的に扱うのは系が巨大なために困難であり、近年計算機の進歩と方法論の発達により扱えるようになったばかりである。しかし、遷移金属を扱うためには非常に高度な電子状態計算を必要とするため、定量的には十分な精度で扱えないのが現状である。
一方、溶液内における電子状態理論において、積分方程式と分子軌道法を組み合わせた
RISM-SCF 法により微視的な溶媒和構造を再現し、物性や反応をよく再現できることが過去の研究により知られている。本研究では、この積分方程式と分子軌道法を組み合わせた方法を上で述べたような酵素などの大きな系に適用し、定量的に十分な精度で求めることのできる新しい方法を開発することを目的とする。今回の発表ではその第一歩として、酵素の活性部位をかたどったモデル化合物の溶液中での物性が、現在の RISM-SCF 法でどれだけ再現できるのか示す。具体的な物質として生体内の電子移動反応に重要な Iron-Sulfur
Protein の一つである[2Fe2S]ferredoxin のモデル化合物（図 1、R=phenyl, o -xylyl などがある）を扱う。この物質の Heisenberg coupling constant、溶液中での励起スペクトル及び
還元ポテンシャル（酸化状態 Fe(III)－Fe(III) から還元状態 Fe(III)－Fe(II) に変化）などの物性を RISM-SCF 法より求め、実験と比較する。
計算は R をメチル基で置き換えたもので行った。用いた溶媒は DMSO である。電子状態
計算は CASSCF、MRMP などを用い、基底は Fe, S には SBKJC の ECP basis に polarizationとさらに S には diffuse を加えたものを、それ以外には 6-31G を用いた。
計算結果の一部を以下に示す。溶液中において酸化状態で構造を最適化した後、励起ス
ペクトルを求めると表１のようになった。真空中と溶媒状態ではスペクトルに変化は見られなかったが、CASSCF に動的電子相関を加えた MRMP では低波長側にシフトすることが分かる。また、それぞれの電子状態を調べたところ、これらのスペクトルが中央の S から Fe への Charge Transfer によるものであることが分かった。この MRMP のスペクトルの値は、実験で得られた 493nm のブロードなピークとほぼ一致する。
また、還元ポテンシャルの値は表２のようになった。なお、CASSCF と MRMP 計算において ∆µ は変化しないものとし、∆E redox = ( E solute (red) + ∆µ (red)) − ( E solute (ox ) + ∆µ (ox ))= ∆E solute + ∆(∆µ )である。これを見ると溶質自身は還元状態（電荷－3）より酸化状態（電荷－2）の方が安定しているが、逆に溶媒和エネルギーは還元状態の方が電荷が大きいため安定しており、還元ポテンシャル全体としては還元状態の方が安定となっていることが分かる。この傾向
は実験と一致しているが、定量的な値（R＝phenyl で 3.50eV、R= o -xylyl で 3.25eV）からはややずれている。これは、実際の物質はベンゼン環がついているために相対的により還元状態の方が安定するためだと思われる。
その他の結果、及び議論は当日発表する予定である。

R

S

S
Fe

R

S

S

R

S

R

Fe
S

図１：モデル化合物
表１：励起スペクトル（単位 cm-1/nm）

Ag→B1uAg→B2uAg→B3uCASSCF/gas
13544/738
13499/741
11548/866

CASSCF/DMSO
13529/739
13568/737
11564/865

表２：還元ポテンシャル（単位 eV）

CASSCF⊿Esolute-4.82⊿(⊿μ(HNC))
6.88
Eredox
2.06

MRMP
-5.13
6.88
1.75

MRMP/DMSO
19622/510
19962/501
16897/592

4Pp016
ガスハイドレートゲスト分子の分子動力学シミュレーション(金沢大理)○松本大輔、井田朋智、加藤信彦、水野元博、遠藤一央【序】ガスハイドレートは水素結合による氷状構造を形成する水分子からなる包接化合物で極地の永久凍土域や水深 500m 以上の海底に広く分布している[1]。自然界では空のガスハイドレートは見つかっておらず、第二の構成要素、ゲスト分子を含むことでその構造が安定化される。ゲスト分
子の大きさにより二種類のハイドレート(Structure I ; sI, Structure II ; sII)構造が見つかっている。sI ハイドレートはゲスト分子が比較的小さい時－例えば二酸化炭素やメタンのような分子－に形成する。より大きなゲスト分子－例えばプロパンのような分子－ではゲスト分子を閉じ込めるためにより大きな空隙を持つ sII 構造が必要とされる。両構造とも単純立方の結晶格子を持っていて、二つの異なるサイズの cage を含んでいる。図 1 に sI のユニットセルを示す。sI のユニ
ットセルには水分子 20 個からなる small cage(12-hedron)二つと水分子 24 個からなる large
cage(14-hedron)六つが含まれている[2]。
ガスハイドレートは体積にして約 180 倍もの天然ガスを閉じ込めることができるので将来のエネルギー源及び液化天然ガスにかわる新しい天然ガスの貯蔵・輸送手段としても期待されている。さらに温室効果ガスの削減方法に
二酸化炭素ハイドレートの生成が含まれており、二酸化炭
素ハイドレートがどのように生成し、成長するかという研究も行われている[3]。このような研究においてはハイド
レートの動的性質の理解が重要となる。特に疎水性のゲス
ト分子がなぜハイドレート構造を安定化させているのかという疑問に答えるため、ハイドレート中でゲスト分子がどのような動的挙動を示すかという問題は興味深い。
これまでハイドレートの動的性質を研究するために分子動力学法が有力な方法として用いられ興味深い結果が得られている[4,5]。しかしこれらの分子動力学シミュレー
ションでは分子間相互作用を表すために古典的な有効ポテンシャルモデルが用いられており、経験的に分子間ポテ
ンシャルのパラメータを与える必要がある。また、このポ
テンシャルは二体相互作用のみ考慮しているため、多体効

small cagelarge cage

図 1 sI のユニットセルと sI を構成する
二種類の cage果を正確に表せない。我々は経験ポテンシャルを用いることなく、また多体相互作用を考慮した分子動力学シミュレーションを行うために原子間(分子間)ポテンシャルを分子軌道計算によって
求める。

【方法】
sI ガスハイドレートのクラスターモデルを考え、ゲスト分子の種類を変えて NVT 一定の分子
動力学シミュレーションを行い、ゲスト分子の運動の違い、small cage と large cage での運動の違い、温度によるゲスト分子の運動の違いを調べる。数値積分には予測子－修正子法を用い、温度を一定に保つために Nose の方法を用いる。分子(原子)間相互作用の計算には半経験的分子軌道計算を用いる。
【文献】[1] E. D. Jr. Sloan, "Clathrate Hydrates of Natural Gases", 2nd ed. Marcel Dekker, Inc.: New
York, 1998.
[2] R. K. McMullan, G. A. Jeffrey, J. Chem. Phys. 42 (1965) 2725.
[3] R. Radhakrishnan, B. L. Trout, J. Chem. Phys. 117 (2002) 1786.
[4] C. Moon, P. C. Taylor, P. M. Rodger, J. Am. Chem. Soc. 125 (2003) 4706
[5] S. Horikawa, H. Itoh, J. Tabata, K. Kawamura, T. Hondoh, J. Phys. Chem. B
101 (1997) 6290

4Pp017
高次元アルゴリズムによる抗不安薬の分子動力学的研究(ＡＴＲ適応研・ＮＴＴ物性基礎研)*○大田原一成、下川信祐、寺前裕之【 はじめに 】
抗不安薬として良く知られたベンゾジアゼピン系誘導体は、ベンゼン環並びに窒素を含む７員環のジシクロ環骨格を持ち、チエノジアゼピン系誘導体は、ベンゼン環がチオフェン環に置き換わった骨格を持つ。この骨格に多様な置換基が付くことにより、多種の薬剤が合成されており、また置換基の違いにより薬効の違いがあることからドラッグデザインの対象として興味深いと考えられる。本研究では、非経験的な計算手法を用いたドラッグデザインを目指して、ベンゾジアゼピン(BZP)系およびチエノジ
アゼピン(TZP)系の抗不安薬について、分子動力学法として高次元アルゴリズムを用いた計算を行った。Ab initio分子軌道法の構造最適化手法に、高次元アルゴリズムを用いることで、コンフォーマーを探索できることは以前に報告した。本報告では、国内で市販されているベンゾジアゼピン系およびチエノジアゼピン系の抗不安薬のほぼ全てにあたる17種類を計算対象とし、薬効の目安である力価と電子状態との関係について考察した。
【 計算方法 】
高次元アルゴリズムによるダイナミクスは、
次の手続きで行った。系のエネルギーは式(1)のハミルトニアンで与えた。ここで pi , p j は運動量、 aij は運動に混合性を導入する行列、 ri は分子内の i 番目の原子の位置座標、V (ri ) は力学系

1
∑ aij pi p j + V (ri ),
2
∂H
r&i =
∑ aij p j ,∂pi j∂H
p& i = −
.
∂ri
H=

(1)
(2)

(3)
のポテンシャルエネルギーで、ab initio 分子軌
道計算で求めた分子のエネルギーを与える。こ
こで系の全エネルギーを保存しながら、式(2),(3)で与えられる運動方程式に従い、ヴェルレ法を用いて各原子の次ステップの位置座標を求めた。
【 結果と考察 】
初期座標を与えた後、ab initio分子軌道法(3-21G基底)を用いて、分子ダイナミクス計算を1000回～3000回行った。最安定構造を探索するため、ダイナミクス100回毎に構造を取り出し、その構造について準Newton法を用いて局所安定構造を得た。図１に例としてクロチアゼパム分子についての繰り返し計算における全エネルギ
ーをプロットした。尚、17種の抗不安薬についての非経験的分子軌道法によるダイナミクス計
算は、非常に大規模であり、PCクラスターを用いて並列計算を行うことで初めて実現された。
PCクラスターは、CPUにPentium4(2GHz)16CPU、 図１．繰り返し計算での全エネルギー

各々のCPU基板にメモリーを1.5GB配し、OSはLinux Kernel 2.4を用いた。
計算対象としたのは、実際に市販されているBZP系およびTZP系の抗不安薬のほぼ全てにあたる17種類で、アルプラゾラム、ブロマゼパム、クロラゼプ酸２カリウム、
クロチアゼパム、クロロジアゼポキジド、クロキサゾラム、ジアゼパム、エチルロフ
ラゼペート、エチゾラム、フルジアゼパム、フルタゾラム、フルトプラゼパム、ロラ
ゼパム、メダゼパム、メキサゾラム、オキサゾラム、プラゼパムの各分子である。
計算したBZP系, TZP系抗不安薬分子の電子状態としてHOMO, LUMO付近のエネル
ギーレベルと、薬効の指標である抗不安薬の力価との関係を調べた結果、HOMO－1のエネルギーレベルが力価との相関が見られた。図２に各抗不安薬分子のHOMO－1
エネルギーと力価の関係をプロットした。HOMO－1 エネルギーが低くなるにつれて力価が強くなる傾向がわかる。現在、HOMO－1 の性質を調べているが、この様な傾向が見出されたことで、非経験的計算手法によって求める電子状態からの抗不安薬デザインの可能性がでてきたと思われる。
1.4
Alprazolam
Fludiazepam
1.2
Lorazepam

Relative Strength1

Mexazolam
Ethyl loflazepate
Flutprazepam

0.8

Etizolam
0.6
Cloxazolam
0.4

0.2

0
-0.37

-0.36

Clotiazepam

Bromazepam
Diazepam
Flutazolam
Prazepam
Dipotassium
clorazepate Oxazolam

-0.35

-0.34

-0.33

Medazepam-0.32

Chlordiazepoxide

-0.31

-0.3

Orbital Energy of HOMO-1 [a.u.]

図２．BZP系およびTZP系抗不安薬の分子軌道エネルギーと力価の相関関係
謝辞
本研究の一部は、通信・放送機構の研究委託により実施したものである。
参考文献
1) 大田原一成, 下川信祐, 寺前裕之,“高次元アルゴリズムによるコンフォーマー探索”, 分子
構造総合討論会2002, 1P130.
2) 新上和正, “高次元アルゴリズム －最適化問題を解く１つの方法－”, 日本ファジー学会誌,
11(3), 382(1999).
3) 寺前裕之, 大田原一成, “高次元アルゴリズムによる分子構造の最適化”, 分子構造総合討論会2001, 327(2001)
4) 大田原一成, 下川信祐, 寺前裕之, “分子構造最適化の並列処理”, 日本コンピュータ化学会
2002春季年会, 2O09(2002).
5) 融道男, “向精神薬マニュアル第２版”, pp.188-198, 医学書院(東京), 2001.

4Pp018

DFT 計算によるトロポロン
二量体における水素移動（京大院理）○今城文雄【序】結晶状態でトロポロンは分子内と分子4'間の二重水素結合による右図のような二量体を形成している（空間群 P21/c, Z=4）1。固体
NMR の測定 2-4 から固体状態での分子運動として、分子の自己拡散、水素移動、180°反転

O2'

3'
2'

OH

O1

O

7
1

6
5

5'
6'

(π-flip) の存在が認められているが、個々の運

1'
7'

O
O1'

H O

2

O2

3

4

動の詳細な情報については不明である。全体として報告されている活性化エネルギーは 24-6
kcal/mol であって、類似した構造をもつ化合物の自己拡散に要するエネルギー（>40 kcal/mol）に比べてかなり小さい。そこで、基礎的データとしてトロポロン二量体自体の水素移動について、その機構とエネルギーを DFT 法により求めたところ、類似構造のカルボン酸二量体とは著しく異なった挙動を示すことがわかった。計算結果を基にして固体状態でのトロポロンの分子運動について新たな機構を提案する。
【方法】分子構造とエネルギーは Gaussian 98 を用いて DFT 法（B3LYP/6-311++G**）により求めた。固体状態における分子間ポテンシャルエネルギーとしては MM2 分子力学計算における van der Waals エネルギーと電荷エネルギーのみを考慮した。
【結果と考察】1. 水素移動.

基底状態 (GS) と遷移状態 (TS) のエネルギーを表 1 に示した。

表 1. トロポロン二量体のエネルギー (E), 零点エネルギー (ZPE),
相対エネルギー (∆E, ∆(E+ZPE))
構造の
対称性
Ci
Cs
D2
C1
D2h
C2v

状態
GS
GS
TS
TS
TS
TS

E
(hartree)-841.81140-841.78846-841.79619-841.78833-841.78620-841.77662∆E(kcal/mol)
0.0
14.4
9.5
14.5
15.8
21.8

ZPE
(hartree)0.22946
0.227550.226170.22195
0.22734
0.21983∆(E+ZPE)(kcal/mol)
0.0
13.2
4.8
12.4
9.8
20.5

二量体の最安定構造の対称性は X 線結晶解析と同じく極めて平面に近い Ci であった。水素移
動の遷移状態はカルボン酸二量体で多く見られる対称性の高い D2h ではなく、二つの芳香環が直交に近い (83°) D2 であった。ただし、完全に直交した C2v は著しく高エネルギーな状態である。水素移動の遷移状態はこの D2 構造以外に一つの分子間水素結合を維持したまま芳香
環がねじれ、二つの芳香環が電荷分離した形の構造が存在することがわかった。ポテンシャ

ル上での遷移状態は C1 で1あるが、実際には二つの芳D2h香環が直交した Cs である。
この水素移動機構では片11一方の芳香環のみが π-flip
することになる。

Ci

1

2

Ci

2

2

2

D2

トロポロンの二量化エ1

ネルギーは BSSE を考慮す
ると 5.3 kcal/mol となり、2カルボン酸に比べ 1/3 程度

C1(~Cs)の大きさしかなかった。こ

2

Ci

1

図 1. 水素移動の機構

れはトロポロンが単量体でも分子内水素結合により安定化することに起因している。
2. 固体中での分子運動. X 線結晶解析と固体 NMR から表 2 に示

A:B

した内の分子運動 1-4 の存在が提

[B-B]

[B-B]

案されている。結晶構造の要請か
ら非等価なサイトへの自己拡散
は水素移動か π-flip が必要とな[A-A]る。今回の計算結果を基にすると、水素移動を伴うが二量体の片一方だけが π-flip する分子運動 5 と図 2. 結晶内での分子配置

上下に芳香環が重なった二分子
A:B（図 2 参照）が擬二回軸のま表 2. 固体状態でのトロポロンの分子運動
a

自己拡散
構造変化
π-flip 水素移動
1. 等価サイト
××[A-A]→[A-A]
2. 等価サイト○○[A-A]→[A’-A’]
3. 非等価サイト○
×[A-A]→[B-B]
4. 非等価サイト
×○[A-A]→[B’-B’]
b
5. なし
○
[A-A]→[A-A’]○c
6. （置換）×
A:B→B:A○aA,B は非等価なサイトを表し、(’) は原子の交換を表す.
b
片側のみ. c 擬 2 回軸まわりの回転.

わりに 180°回転する分子運動 6の存在が考えられる。分子運動 5と 6 は実質的には分子の自己拡散を伴わないので、分子運動 1-4よりも有利であると考えられる。
実際、分子運動 5 に必要なエネル
ギーは二量体に関しては実験値
よりも十分に小さい。また、分子運動 6 について、剛体モデルを用いて A:B 間の相互回転などを考慮した分子間ポテンシャルを計算したところ、その障壁は二量化エネルギーを含めても 61 kcal/mol 程度であった。分子間のより精密な最適化を行えば容易に実験値に近い値になると考えられる。
1. H. Shimanouchi and Y. Sasada, Acta Crystallogr., B29, 81 (1973).
2. N. M. Szeverenyi, A. Bax, and G. E. Maciel, J. Am. Chem. Soc., 105, 2579 (1983).
3. A. Detken, H. Zimmermann, U. Haeberlen, and Z. Luz, J. Magn. Reson., 126, 95 (1997).
4. Z. Gan and R. R. Ernst, J. Chem. Phys., 108, 9444 (1998).

4Pp019
メタノール水素結合錯体の振動緩和（九大総理工）○ 西野裕美 ・ 山本典史 ・ 三好永作
■ 序論：凝縮相で観測される振動エネルギーの散逸過程（振動緩和）では、反応系に不特定多数の溶質・溶媒分子が関与しているので、全ての相互作用を厳密に考慮した解析は実験的にも理論的にも困難な課題である。一方、孤立気相で観測される
水素結合錯体は、特定少数個の溶質・溶媒分子のみで構成されているので、分光実験や理論計算から安定構造を比較的精度良
く特定することが可能となる。このような水素結合錯体は、溶質分子内や溶質・溶媒分子間で起こる振動緩和を分子レベルで調査するための重要な手掛かりと考えられる。本研究では「振動緩和における溶媒和分子の役割を明らかにする」ことを目的
として、メタノール分子および水和メタノールを対象に、メタノール OH 伸縮モードを振動励起した後に起こるエネルギー再分配の時間発展過程を、古典トラジェクトリー解析を用いて調査している。

■ 計算：HF/6-31G(d,p)構造最適化で得られたメタノール分子・水和メタノールの各平衡構造を初期構造として、OH 伸縮
モード（図１d）の n 倍音（n=1-5）を振動励起した後のエネルギー再分配の時間発展過程について、古典トラジェクトリー
計算を用いて追跡した。今回用いた古典トラジェクトリー計算では、各時間ステップで HF/6-31G(d,p)計算を実行すること
で各原子核に働く力を正確に求めている。得られた計算結果は、基準振動モード（図１）を用いて解析を行った。

■ 結果：OHS 伸縮（OHS）モードの５倍音を励起して得られた結果を示す。図３は 128 fs 経過後の OH 変角（OHB）モード、CH 伸縮（CHS）モードおよび OHS モードを含む配位空間上の代表点の軌跡を示す。図４は各振動モードに対応する運動量成分の時間変化を示す。ここで、図中のグラフはそれぞれ、OHB モード（上段・青）、CHS モード（中段・緑）および OHB モード（下段・赤）の運動量成分である。図５は 128 fs 経過後の位相空間上の代表点の軌跡を示す。

（ａ）OHB モード（ｂ）CHB モード（ｃ）CHS モード（ｄ）OHS モード

図１．水和メタノールの基準振動モード

0.0

-0.2

-0.4
-0.4

-0.2

0.0

0.2

0.4

0.6

Q̲OHS(Bohramu**1/2)

0.2

0.0

-0.2

-0.4
-0.4

-0.2

0.0

0.2

0.4

Q̲OHS(Bohramu**1/2)

Q1=OHB, Q2=OHS

Q1=CHS, Q2=OHS

0.6

0.4

Q̲OHB(Bohramu**1/2)

0.2

0.4

Q̲OHB(Bohramu**1/2)

0.4

Q̲OHB(Bohramu**1/2)

Q̲OHB(Bohramu**1/2)

0.4

0.2

0.0

-0.2

-0.4
-0.4

-0.2

0.0

0.2

0.4

Q̲OHS(Bohramu**1/2)

0.0

-0.2

-0.4
-0.4

-0.2

0.2

0.4

Q1=CHS, Q2=OHS（ｂ）水和メタノール

図２．ポテンシャルエネルギー曲面

0.0

Q̲OHS(Bohramu**1/2)

Q1=OHB, Q2=OHS（ａ）メタノール分子

0.6

0.2

0.6

ａ．メタノール分子の場合：運動量の時間変化（図４a）をみると、OHS モードが励起した後、すぐに CHS モードへとエネ
ルギーが流れ込むことが分かる。しかし、このとき分配されるエネルギーはほんの僅かであり、計算過程で全ての振動モード
間の相互作用を正確に取り込んでいるにも係わらず、実際には CHS モード以外の振動モードにはエネルギーがほとんど分配されない。更に、256 fs 程度時間が経過すると、与えられたほとんどすべてのエネルギーが最初に励起された OHS モードに再帰することが明らかとなった。このように、メタノール分子の場合、OHS モードの５倍音を励起した後の振動ダイナミ
クスでは、体系の状態が周期的に初期状態へ再帰する現象（再帰現象）が観測された。また、配位空間（図３a）上の代表点についても規則的な軌跡を描いて、周期的に初期状態 Qi=0 へ再帰している。

ｂ．水和メタノールの場合：OHS モードを励起した後、メタノール分子と同様に、CHS モードが励起される（図 4b）。励起された CHS モードは、与えられたエネルギーがその後しばらくはほとんど増大しないのに対して、64 fs 程度経過すると、OHB モードが急激に大きなエネルギーを持つようになる。このとき、メタノール分子で観測された再帰現象は現れず、時間が経つにしたがって様々な振動モードが次々と励起される。その結果、メタノール分子では、与えられたエネルギーが少数の振動モードのみに局在していたのに対して、水和メタノールでは、与えられたエネルギーが多数の振動モードに非局在化して．．．
いる。位相空間上の代表点の軌跡（図５b）をみると、OHS モードと CHS モードは運動量を徐々に変化させながら振動運動．．．
しているのに対して、OHB モードは、一定時間経過した後で、運動量を急激に変化させながら振動運動するようになることが分かる。このことから水和メタノールでは、OHB モードが振動緩和に対して重要な役割を果たしていると考えられる。

当日の発表では、他の n 倍音（n=1-4）で得られた結果や今回取り上げた以外の振動モードの挙動も示しながら、様々なメタ

0.4

0.4

0.3

0.3

0.3

0.2
0.1
0.0
-0.1
-0.2

0.2
0.1
0.0
-0.1
-0.2

Q̲CHS(Bohramu**1/2)

0.4

0.3

Q̲OHB(Bohramu**1/2)

0.4

Q̲CHS(Bohramu**1/2)

Q̲OHB(Bohramu**1/2)

ノール水素結合錯体に関する数値計算の結果をもとにして、振動緩和における溶媒和分子の役割について議論する予定である。

0.2
0.1
0.0
-0.1
-0.2

-0.3

-0.3

-0.3

-0.4
-0.4

-0.4
-0.4

-0.4
-0.4 -0.2 0.0 0.2 0.4 0.6
Q̲OHS(Bohramu**1/2)-0.2 0.0
0.2
0.4 0.6
Q̲OHS(Bohramu**1/2)

0.8

Q1=OHB, Q2=OHS-0.2 0.0
0.2
0.4 0.6
Q̲OHS(Bohramu**1/2)

0.8

Q1=CHS, Q2=OHS

0.2
0.1
0.0
-0.1
-0.2
-0.3
-0.4
-0.4 -0.2 0.0 0.2 0.4 0.6
Q̲OHS(Bohramu**1/2)

0.8

Q1=OHB, Q2=OHS（ａ）メタノール分子

0.8

Q1=CHS, Q2=OHS（ｂ）水和メタノール

0.1

0.0

0.0

-0.1
-0.2
0.2
0.1
0.0
-0.1
-0.2
0.4
0.2
0.0

-0.1
-0.2
0.2
0.1
0.0
-0.1
-0.2
0.4
0.2
0.0

-0.2

-0.2

-0.4
0.0

-0.4
0.0

64.0

128.0

192.0

256.0

0.4

0.3

0.3

0.2
0.1
0.0
-0.1
-0.2

64.0

128.0

192.0

t(fs)（ｂ）水和メタノール

図４．運動量の時間変化

256.0

0.2
0.1
0.0
-0.1
-0.2
-0.3

-0.3

t(fs)（ａ）メタノール分子

0.4

P(Bohramu**1/2fs**-1)

0.2

0.1

P(Bohramu**1/2fs**-1)

0.2

P(Bohramu**1/2fs**-1)

P(Bohramu**1/2fs**-1)

図３．配位空間

-0.4
-0.4

-0.2

0.0
0.2
0.4 0.6
Q(Bohramu**1/2)（ａ）メタノール分子

0.8

-0.4
-0.4 -0.2

0.0 0.2 0.4 0.6
Q(Bohramu**1/2)（ｂ）水和メタノール

図５．位相空間0.8

4Pp020ab initio 法による水クラスターイオン[H+(H2O)n (n=2-6)]の核生成過程に関する理論的研究（早大理工 a ・理研 b）○小澤志保 a、山内佑介 a、河井葉子 b、山口悟 b、
岡田芳樹 b、中井浩巳 a【序】結晶成長初期過程における核生成現象は典型的な非平衡過程である。酸性雨発生機構などの自然現象のみならず、半導体生成プロセスや生体反応などに関連していることから重要視され、理論と実験の両面から様々な研究が行われている。しかしながら、核生成進行の挙動についての理解はいまだ不十分である。
均一核生成に対して、コアがイオンであるイオン核生成のプロセスはマイクロクラスターに準安定状態を持ち、反応速度が著しく増大することが知られている。特に水クラスターイオン H+(H2O)n は大気中に広範囲に存在し、様々な分子と反応して酸性雨発生機構などに関わる物質である。理研の岡田らはサイズ選別したクラスターイオンを様々な分子に衝突させる実験を行った。
水クラスターイオン＋水：
水クラスターイオン＋アセトン：
水クラスターイオン＋DMSO：

H+(H2O)n+ D2O

(1)

+

(2)

+

(3)

H (H2O)n+ CH3COCH3
H (H2O)n+ CH3SOCH3

その結果、イオンにモノマーが取り込まれる吸着、クラスターが分解してしまう解離、という 2 つの反応過程が存在することが明らかになり、その反応断面積の値は、衝突エネルギーやクラスターサイズの変化に対して、それぞれ特徴的な変化を示すことがわかった。しかしその違いについての原因は明らかでない。
本研究では各反応における引力的ポテンシャルに基づく衝突断面積を、ab initio 法を用いて理論的に求める。この結果から各反応の特徴やマイクロクラスターについて議論する。
また衝突反応の結果、吸着が起こったクラスターは高い振動励起状態にあると考えられる。結晶成長が安定して進行するためには、この高い振動励起状態から脱励起する必要がある。そこで我々は衝突反応とその後に続く冷却過程を ab initio molecular dynamics (AIMD) シミュレーションによって検討した。

【水クラスターイオン[H+(H2O)n (n=2-6)]核生成過程の解析的衝突断面積】
距離 R、並進エネルギーEcol および衝突パラメータ b が与えられているとき、衝突に対する障壁として遠心障壁を加味した有効ポテンシャル Veff = V(R) + Ecol b2/R2 を考える。この時 Veff は遠心斥力と V(R)の遠距離性引力部分の組み合わせの結果として生じ、R = Rmax で局所的な最大値をもつ。衝突断面積を求めるためには以下の式を陰関数的に解き、Rmax および bmax を求めればよい。


1 &2
b2 
µ ( R ) R max =  E col − V ( R) − E col 2 
≥ 0,
2
R

 R max

(4)


b2 
E
−
V
(
R
)
−
E
= 0.
col
max
 col
2 
Rmax

 b max(5)

これより、衝突断面積 σ は σ= πbmax 2 として得ることができる。
図 1 および図 2 に結果を示す。まず水クラスターイオン[H+(H2O)n]の H3O+を中心とする構造に対して衝突断面積を求めた。n=4 までは中心の H3O+が H2O により遮蔽されるが、衝突断面積は緩やかに増加している。結果として H3O+を中心とする構造での衝突断面積を求めた。この構造は ab initio 計算では
n=2,4,6 に安定構造の存在が予想できた。実験の振動スペクトル解析では n=2,6 に H5O2+を中心とする安定構造の存在が予測されている。ここで H5O2+を中心とする構造の衝突断面積は、H3O+を中心とする構造と比べ n=6 で減少しており、実験の挙動と一致している。結果として n=6 では H5O2+を中心とする構造をとることが断面積という物性から予測できた。

水クラスターイオンにアセトンまたは DMSO を衝突させる衝突断面積は、水モノマーを衝突させた場合と似た挙動を示している。その反応確率は水、アセトン、
DMSO の順に大きくなっており、その傾向はプロトン[Å2]
Cross section

水クラスターイオンの H5O2+を中心とする構造の衝突◇H3O+center○H5O2+center

Calc.

400
300

Exptl.

200
100
1

2

3

4

5

6

7

6

7

Size
Fig. 1. Calculated and experimental crosssection of H+(H2O)n + H2O (n=2-6).
1000

[Å2]

験の反応断面積をある程度再現している。続いて、

500

Cross section

断面積は n の増加に伴って単調増加を示し、衝突実

600

DMSO

900
800

Acetone

700
600
500
1

親和性に比例していることが分かった。

2

3

4

5

Size
Fig. 2. Calculated collision cross section of
H+(H2O)n + X (n=2-6), where X indicates
Acetone or DMSO.

【水クラスターイオン[H+(H2O)n (n=2)]振動緩和のAIMD シミュレーション】

高い振動励起状態にある水クラスターイオンが、不反応分子の衝突により脱励起される過程を考察する。
そのために以下の反応を AIMD シミュレーションによって扱った。
H+(H2O)１ + H2O → H+(H2O)2*
+

*

+

(6)

H (H2O)2 + N2 → H (H2O)2 + N2*

(7)

ここで＊は振動励起状態にある分子を表す。(6)の反応によって振動励起された H+(H2O)2*が生成するが、これは分子振動緩和(IVR)および分子間振動エネルギー移動の 2 種類の過程によって緩和される。そこで衝突後の H+(H2O)2 のトラジェクトリーから IVR の進行が異なる 4 つの状態を取り出し、(7)の反応におけ4 つの状態に対し、N2 を様々な配向から衝突させた。シミュレーションから動の様子を観察した。図 3 は 1 つのト
ラジェクトリーにおける H+(H2O)2 と N2の全エネルギーの時間変化を示して+いる。衝突に伴い H (H2O)2 の全エネ
ルギーが N2 へ移動している様子が見
られる。

25

60
40

20

E N2

20

15

0
-20

10

E H+(H2O)2

-40

5

R

0

-60
0

0.5

1

Time [ps]1.5

2Fig. 3. Time evolution of total energy for
H+(H2O)2 and N2 (E H+(H2O)2, E N2), anddistance between H+(H2O)2 and N2 (R).

Distance [Å2]

H+(H2O)2 と N2 の間でのエネルギー移

Total energy (E kine+E pot) [kcal/mol]

る H+(H2O)2*の初期条件とした。その

4Pp021アセチレン類似ケイ素およびゲルマニウム化学種と遷移金属錯体の相互作用に関する理論的研究(京大院工)○蔵本有紀・澤井則行・藤原雄介・隅本倫徳・
中尾嘉秀・佐藤啓文・榊 茂好【序】Si2H2 および Ge2H2 はアセチレン型の三重結合を持つ構造(1)よりも、H 1,2-移動型（ビニリデン型; 2）や 2H 架橋型(3)が安定であり、最も安定な構造は 3 である（Scheme 1 参照）。これらの相対安定性および異性化の活性障壁は理論的にも良く研究されている 1)。一方、有機金属錯体の化
学では、アセチレンやビニリデン錯体の存在が古くから知られており、それらの構造と結合性の理論的解明も進んでいる。
従って、Si2H2 や Ge2H2 の遷

123.8

127.7
2.182
1.478

1.680
49.2
2.198
72.5

1.475

2.074

1.625

1.744
54.1
2.097

163.1
1.478

移金属錯体の存在も同様に期待され、その構造や結合

1

2

3

4

性、相対安定性がフリーのScheme 1

Si2H2 や Ge2H2 とどう異な

DFT 法による最適構造

るのか興味深い。しかし、Si2H2 の遷移金属錯体の理論的研究はこれまで行われていない。本研究では Pt(PR3)2(Si2H2), RhCl(PMe3)2(Si2H2), Cp2Zr(Si2H2)の理論的研究を行い、構造と相対安定性、結合性を検討した。
【計算方法】
構造最適化は DFT 法を用いて行った。交換相関項は B3LYP 汎関数で近似した。Pt, Rh, Zr, Si, Ge,
P, Cl の内殻電子は有効内殻ポテンシャル(ECP)で置き換え、それらの価電子には split valence 型基底関数を用いた。その他については 6-31G*を用いた。結合性、相対安定性の検討には、より良い基底関数系を使用した。Pt, Rh, Zr の内殻電子は構造最適化計算と同様に ECP で置き換えたが、Ptの原子価電子は(541/541/111/1)で表し、他の金属も同様の基底関数で表した。Si2H2, Ge2H2 には6-311G(d,p), 6-311G(2d,2p)を使用した。
【結果と考察】表１．Si2H2 の相対エネルギー(kcal/mol)

最初に、フリーの Si2H2 のSi2H2 の基底関数

1

2

3

4

LanL2DZ

13.5

5.8

0.0

6.8

6-311G(d,p)

18.6

10.5

0.0

9.8

6-311G(2d,2p)

20.2

11.2

0.0

10.5

MP4(DQ)

6-311G(2d,2p)

23.4

12.7

0.0

12.7

MP4(SDQ)

6-311G(2d,2p)

22.0

12.6

0.0

12.3

MP4(SDTQ)

6-311G(2d,2p)

18.8

13.8

0.0

10.8

CCSD

6-311G(2d,2p)

21.4

11.6

0.0

12.2

LanL2DZ

8.7

5.0

0.0

6.0

6-311G(d,p)

16.6

12.4

0.0

10.6

6-311G(2d,2p)

18.4

12.3

0.0

10.9

各種の構造の相対安定性が、どの計算レベルで正しく再現されるか、検討した。構造は DFT 法で最適化されたものを採用した。表 1 に示すよ
うに、相対安定性に対するSi の基底関数効果は大きく、
LanL2DZ ではエネルギー差

DFTが過小評価され、6-311G(d,p)
および 6-311G(2d,2p)では同

CCSD(T)

程度の相対安定性を与える。しかし、1 の相対安定性は 6-311G(d,p)と 6-311G(2d,2p)でなお 2 kcal/molの差が見られる。また、3 に対する 1 の相対安定性は MP4(SDQ)と MP4(SDTQ)の間、CCSD とCCSD(T)の間でかなりの差を示し、3 電子励起の考慮が必要であることを示している。今後、より詳細な検討を行う予定である。DFT(B3LYP)法は MP4(SDQ), CCSD と MP4(SDTQ), CCSD(T)の中間であり、相対安定性をある程度の精度で再現している。
Figure 1 に Pt(PMe3)2(Si2H2)の最適構造を示した。1-Pt では
Si-Si 距離は 0.06 Å 長くなっている。2-Pta では Si=Si 結合が

108.9

2.370

Pt 2.370

2.385

2.385

1.485 Si

side-on 配位しているが、2-Ptb

P

P

2.398

2.400

2.446

1-Pt

では end-on 配位している。

P

110.6

2.425 Pt
1.486

Si Si-Si-H = 135.9 °

2.134

P

P

Si

2.285

Si

長くなっているが、2-Ptb では

2-Pta

109.7

Pt 2.385

Si

Si-Si 距離が 0.17 Å 長く、4-Pt

P

P
2.385
2.438

ほとんど変化がない。3-Pt では

2.438

2-Ptb
P

109.9

P

2.367

Pt

2.373

2.316

2.366

1.485

Si
1.685

2.577

Si

2.203

Si

Si-Si-H = 145.4 °

1.728

1.643

3-Pt

では約 0.1 Å 長くなっている。

2.209
Si
2.186 117.5
1.480 Si
Si-Si-Pt = 153.9°

Si-Si-H = 121.4 °2-Pta では Si-Si 距離が約 0.1 Å

114.1 P

2.350 Pt 2.350

4-Pt

Figure 1. Pt(PMe3)2(Si2H2)の各種の最適構造

これらの変化はアセチレンやエチレン錯体同様に

表 2．Pt(PMe3)2(Si2H2)および RhCl(PMe3)2(Si2H2)の異性体の相対エネ

Pt(0)から、Si2H2 のπ*軌道

ルギー(kcal/mol)と Si2H2 のπ*軌道のエネルギー準位ε(π*) (eV)への電荷移動が配位結合に含まれていることを示唆している。RhCl(PMe3)(Si2H2)についても、同様の構造が得ら
れている。

1-M

2-Ma

2-Mb

3-M

4-M

Pt

6.3

0.0

12.1

0.9

4.1

Rh

18.1

0.0

17.5

6.1

16.3

ε(π*)

-3.4

-2.0

-2.1

-3.7

DFT 法で求めた Pt, Rh 錯体の相対安定性を表 2 に示した。Pt, Rh 錯体になると、相対安定性に対する基底関数効果はフリーの Si2H2 に比べて小さい。また、DFT 法と MP4(SDQ)法による相違も小さくなる。最も注目すべき結果の一つは、Pt(0)錯体および Rh(I)錯体ではビニリデン型構造が最も安定になることである。また、アセチレン型構造はフリーの Si2H2 と同様に 2H 架橋型構造よりも不安定であるが、そのエネルギー差は小さくなっている。
NBO 解析より、Pt(0)錯体では Si2H2 の電荷密度が増加し、Pt の d 軌道密度は減少しており、Pt(0)の d 軌道から Si2H2 への電荷移動が起きている。Rh(I)錯体では Si2H2 部分の電荷密度が減少し、Rh の d 軌道電子密度は増加しており、π 逆供与の寄与が小さいことが示唆される。フリーの Si2H2で最も安定な 3 では 1 や 2 に比べπ*軌道が高く、Lewis 塩基性の高い Pt(0)や Zr(II)錯体などでは相対安定性が低下したが、Rh(I)錯体では逆供与の寄与が小さく、1-Rh はかなり不安定なままであると考えられる。
【参考文献】
1) S. Nagase, K. Kobayashi, and N. Takagi, J. Organomet. Chem., 611, 264 (2000).

4Pp022[Re2Cl8]2-分子系の基底・励起状態に関する理論的研究(東大院工 a、JST PRESTb) ○羽根慎太郎 a、中嶋隆人 a,b、平尾公彦 a【序】1965年、[Re2Cl8]2-分子系においてはじめて、δ結合を含む金属間の四重結合が見つかった1。
この四重結合は1つの σ結合(Reの5dz2軌道から生成)、2つの π結合(Reの5dxz軌道から1つ、5dyz軌道から1つ生成)、および1つの δ結合(Reの5dxy軌道から生成)から成り、Reの5dx2-y2軌道はClとの結合に使われている。図1に[Re2Cl8]2-分子の構造を、図2に[Re2Cl8]2-分子系の電子スペクトルを示す2。

1.8eV

図 1. [Re2Cl8]2-の構造

図 2. [(η-C4H9N)2Re2Cl8]結晶の電子スペクトル

これまで、電子スペクトルの同定を行うためにさまざまな実験や計算が試みられてきた。しかし、その同定結果は各報告において異なっている。特に、1.8eVに観測される最初の大きなピークの同定については、実験と計算の間に顕著な違いがある。このピークは実験では δ- δ* (1A1g-1A2u)
許容遷移と同定されるが、計算では0.87eV3または3.20-3.384eV4,5と算出されてきたのである。これまで用いられてきた計算方法は拡張Huckel法、X α法、ECPを用いたab initio計算などであった。
All-electronのab initio計算は最近になってようやく1つ報告された。この計算では、δ- δ* ( 1A1g-1A2u)
遷移は2.03eV6と算出された。そこで本研究では、[Re2Cl8]2-分子系に対してより高精度なall-electronのab initio計算を行い、電子スペクトルの同定を試みる。また、同様の計算を類似の系である[Mo2Cl8]4-分子系と[Tc2Cl8]3-分子系についても試みる。
【計算方法】構造は X 線回折による実験値を用いた。重原子を含む系なので、相対論効果を 3 次の Douglas-Kroll(DK3)法により取り入れる。また、電子相関効果を CASPT2 法により考慮する。
基底関数は本研究室で開発された相対論的基底関数を用いた。その exponent を用い、分子計算に適用するため、すべての原子について比率 2.5 の even-tempered によって 0.01 を下らないところま
で exponent を増やした。これと DK3 法により最適化した ANO タイプの縮約係数を用いた。得られた基底関数は、Re: (29s26p15d12f6g)/[8s7p5d3f1g]、Cl: (18s13p6d)/[5s4p2d]となる。CASSCF 計算の active space には、14 軌道に 20 電子を配置(20e14o)する。この active space には Re-Re の結合に使われる 5dz2、5dxz、5dyz、5dxy 軌道すべてと、Cl の 3p 軌道の一部が含まれる。CASPT2 計算のcore 軌道は、Re は 3d 軌道まで、Cl は 1s 軌道までとした。

【結果と考察】まず Hartree-Fock 計算を行ったが、四重結合を記述することはできなかった。また、B3LYP による構造最適化計算を行ったところ、Re-Re の距離が実験値の 2.222Åに対し、2.521Åと大きめに算出された 7。以上の 2 点から、四重結合の記述には多配置計算が必要であると示唆される。これは、過去の研究でもいわれてきたことである。これに対して、CASSCF、CASPT2
計算では四重結合を記述することができた。基底状態に大きく寄与する配置は、Hartree-Fock 配置[…( σ) 2( π) 4 (δ) 2 ( δ* )0…](66%)、HOMO-LUMO の 2 電子励起配置[…( σ)2( π) 4 ( δ) 0 (δ* ) 2…](19%)
であった。表 1 に[Re2Cl8]2-の spin-free レベルにおける CASSCF、CASPT2 計算による d-d 遷移エ
ネルギーの一部を過去の研究と併せて示す。
表 1. [Re2Cl8]2-の spin-free レベルにおける d-d 遷移エネルギー(eV)
Exp.
1

Hay4

Blaudeau5

Gagliardi6

本研究

本研究

Xα

GVB-CIMCSCFCASPT2CASSCFCASPT2

0.87

3.20

3.384

2.03

3.03

1.69

0.45

0.75

0.435

0.40

0.32

3.90

7.48

5.69

3.83

3.54

3.84

3.64

2.98

1

A1g→ A2u(δ-δ*)

1

Bursten3

1.8

A1g→3A2u(δ-δ*)

1

A1g→1A2u(π-π*)

1

4.9

A1g→3A2u(π-π*)

4.56

まず δ- δ* (1 A1g-1A2u)遷移エネルギーに注目する。この遷移はこれまで、Bursten の計算では 0.87eV、Hay や Blaudeau の計算では 3.2-3.38eV と、いずれも実験値との違いは 1eV 以上であった。しかし、本研究においては実験値との違いは 0.1eV となる。これは最近になって報告された Gagliardi による all-electron の計算と同じ同定結果となる。しかし、本研究の結果と Gagliardi の結果には違いもいくつかある。その 1 つである π- π* ( 1A1g→1A2u)遷移エネルギーは、本研究では 3.83eV となるのに対し、Gagliardi の研究では 4.56eV である。この原因は、active space の違いにある可能性が高い。
図 2 に示すスペクトル全体の spin-orbit coupling まで含めた同定結果と詳しい考察、および[Mo2Cl8]4-分子系と[Tc2Cl8]3-分子系の結果については当日発表する。
【参考文献】(1)F. A. Cotton

Inorg. Chem. 4, 334 (1965)(2)W. C. Trogler, H. B. Gray

Acc. Chem. Res. 11, 232 (1978)(3)B. E. Bursten, F. A. Cotton et al(4)P. J. Hay

J. Am. Chem. Soc. 105, 3082 (1983)

J. Am. Chem. Soc. 104, 7007 (1982)(5)J-P. Blaudeau, R. B. Ross et al(6)L. Gagliardi, B. O. Roos

J. Phys. Chem. 98,7123 (1994)

Inorg. Chem. 42, 1599 (2003)(7)F. A. Cotton, B. A. Frenz et al

J. Am. Chem. Soc. 98, 2768 (1976)

4Pp023
ジアジリンの熱分解反応に関する研究
分子研１、筑波大２、Prince Edward Island 大３○崔隆基１、永瀬茂１、赤阪健２、Michael T. H. Liu３
１． 序
ジアジリンは 1960 年 Paulsen らによって報告されて以来、取り扱いの容易さや化合物の安定性などの理由によりカルベン生成の先駆物質として有機合成に幅広く用いられている。ジアジリンは光や熱によりカルベンと N２へ分解またはジ
アゾアルカンへの異性化がおきることが知られている（図１）。熱分解反応メカ
ニズムは、非常に短寿命の反応中間体の検出の難
R
R'
+ N2
C
しさなどの理由により、過去に様々な研究が行わ R R'
C
れてきたにもかかわらずいまだに解明されていR'ない。本研究は計算と実験とのインタープレイに N N
CNN+よりジアジリンの熱分解反応機構を明らかにすRるのを目標とする。
図１
２． 方法
本 研 究 で はphenylchlorodiazirine(1),phenyl-n-butyldiazirine(2),2-admantane-2,3’-[3H]diazirine(3) の 熱 分 解 反 応 を 研 究 し た 。 計 算 方 法 は
B3LYP/6-31G(d)法を用いた。またジアゾアルカンへの異性化反応の遷移状態は多参照波動関数を用いないとうまく記述できないので MRMP(8,8)/CASSCF(8,8)
計算も並行して行った。
３． 結果
図２は 1 の熱分解反応のポテン
シャルを表したものである。斜36.2
体の数字は MRMP(8,8)法によって得られたエネルギーで、その
30.9
18.2
Cl-C-Ph + N2
30.1
他は B3LYP 法で得られた値であ
4.9
4.9
4.0
0.0
る。MRMP レベルでカルベンへ
Cl-C-Ph + N2-7.3の分解反応がジアゾアルカンへの異性化より 6.1kcal/mol 低いのが分かる。実際実験を行ったと
ころカルベンが唯一の生成物と
図２
して検出され、計算結果と良い一致を示している。2 と 3 の結果及び考察は当日報告する予定である。

4Pp024 シラシクロブテンの開環反応に関する理論的研究(九大先導研・倉敷芸科大) ○安永真・塩田淑仁・
蒲池‫ݗ‬志・仲章伸・石川満夫・吉澤一成【序】シラシクロブテンは熱反応で開環反応を֙こし、炭素−ケイ素二重結合を含む不安定中間体を生じる。その際、４員環の開裂のしかたによりシラシクロブテン(1)は 2 種໸の中間体(2、3)を生じる（スキーム１）。しかし、メタノール封管中の実験（160℃）では 2 に由来する生成物のみが報告されている。本研究では、密度汎関数理論を用いてシラシクロブテンの開環反応について考察し、開環反応の選択性を遷移状態の考察により明らかにした。さらに、中間体の炭素−ケイ素二重結合へのメタノール付加反応もӂ析した。

スキーム１
tBu

H

tBu

H

tBu
Me 3SiOSiMe 3Si
Si

OSiMe 3tBu

Me 3Si
SiMe 3

Me 3Si

1

tBu

H

tBu

Si
OSiMe 3

Me 3Si
SiMe 32

3【‫ٽ‬算方法】‫ٽ‬算方法はＢ３ＬＹＰ法、基底関数は 6-31Ｇ*基底を用いた。反応物、遷移状態、中間生成物及び生成物の構造最適化を実行した。安定構造、および遷移状態の構造を最適化した後に振動ӂ析を行い、安定構造、遷移状態であることを確認した。‫ٽ‬算プログラムにはＧＡＵＳＳＩＡＮ98
プログラムを用いた。

【結果及び考察】1 の開環反応は立体障壁のため反旋で進行するのは不可能であるため、同旋で進行する経路のみを考えた。図１に示すエネルギーダイアグラムのように、シラシクロブテンの開環反応の活性化エネルギーは tBu 基がトランス位になる反応経路では 21.2kcal/mol であるのに対して、

tBu 基がシス位になる反応経路では 36.8kcal/mol となり、前者のほうが 15.6kcal/mol 有利である
ことが明らかとなった。また、2 と 3 のポテンシャルエネルギーの差はあまりないことから、1 と中間体が熱的な平衡状態に達しているならば 2 と 3 が同程度存在する可能性がある。しかし、160℃の条件では 3 に対応する生成物は得られていないので、3 にԁべてエネルギー的に有利な 2 のみが発生し、それがメタノールと反応して、上述のような実験結果になると考えられる。いっぽう、実験条件
を不活性ガス中で 190℃の条件で行うと、炭素上の SiMeO 基がケイ素上に転移した副生成物が得られるようになる。立体化学的に SiMeO 基を外側にもつ 2 からはそのような転位は不可能であると考えられるので、この場合には開環反応でエネルギー的に不利な経路を通過して 3 ができていると推測
される。

tBu

H
3
1.7
4

117.8
127.9
7
3.50

Si

6
1.7
4

1.897

tBu

1.491

2

2.016

Me 3SiSi

128.1 123.0
3
3.66

Si
OSiMe 3

5

1
0.0

3
1.3

1.572

104.9

tBu

H

12.3

1.338

96.2

3

13.1

tBu

H

φ(SiCCC) = 68.5

SiMe 33
14.9

tBu
OSiMe 3

Me 3Si

TS1
21.2

3

TS2
36.8

4
1.3

1.496

OSiMe 3tBu

Me 3Si

SiMe 3 φ(SiCCC) = –79.3

1
Me 3Si2

図１開環反応のエネルギーダイアグラム

図２に中間体 2 及び 3 に対するメタノールの付加反応の結果を示す。2 に対するメタノールの付加反応では活性化エネルギー10.1 kcal/mol と小さく、反応熱は 55.0 kcal/mol と大きい。そのため、2 はメタノールの存在下では速やかに反応すると考えられる。3 に関しても 2 の場合と同様の結果がえられ、活性化エネルギーと反応熱は２つの中間体でほとんど差がないということが明らかとなった。つまり、最初の開環反応が律速段階であり、その中間体の選択性が生成物比を決定している。
他のシラシクロブテンでも同様のӂ析を行った。詳細は当日発表する。

(a)

TS3
10.1

(b)
3 + MeOH

2 + MeOH
0.0

TS4
8.9

0.0

H

H

tBu

MeO

MeO
OSiMe3SitBu

Me 3Si
4

–55.0

図2

tBu

Si
OSiMe3

Me 3Si
SiMe 3SiMe 3

4

tBu

H

H

5

5

–52.0

メタノール付加反応のエネルギーダイアグラム【参考文献】
Naka,A.;Ikadai,J.;Motoike,S.;Yoshizawa,K.;Kondo,Y.;Kang,S-Y.;Ishikawa,M.Organometallics20 02,21,2033.
Yoshizawa,K.;Kondo,Y.;Kang,S-Y.;Naka,A.;Ishikawa,M.Organometallics20 02,21,3271.