2Ta01

多配置 SCF 法の新しい可能性（九大院理）中野 晴之

CAS (complete active space)によらない一般の MC-SCF 法の新しい可能性と適用例について発表する．CAS-SCF (Complete Active Space SCF)法は，反応過程の結合の組み換えを正しく記述すること，反応にともなう分子構造の変化に対して安定であること，size-consistent である
ことなどの特長から，化学反応のポテンシャル・エネルギー曲面の記述に際して，非常によく利用されている．しかしながら，しばしば指摘されるように，CI 空間の次元は active 軌道
数に対し急速に増加し容易に計算の限界に達することから，適用範囲が限られてもいる．一方，CAS に属しない一般の MC-SCF 法は，SCF の収束が遅いこと，size-consistency を満たさないことなどから，CAS-SCF 法の蔭に隠れ，これまではあまり使われてはこなかった．しかしながら，近年の SCF の収束法の発展，計算機の速度向上により，SCF の収束の問題は大きく改善された．また，現在では，発表者らの開発した GMC-QDPT などにより，一般の多配置波動関数を出発点として動的電子相関を取り入れ，化学的精度を得ることができるようにもなっており，多配置 SCF 法の新しい可能性がひろがりつつある．本発表では，(P+SD)-SCF
法，あるいは，string-product space (SPS) SCF, occupation specified space (OSS) SCF など，やや特殊な MC-SCF 法を例に，多配置 SCF 法の新しい可能性を考察する．
現在，さまざまな CI 法の計算法が知られている．我々の general MC-SCF プログラムの CI部分では，α-string と β-string の積からなる Slater 行列式を基底として使用している．σ ベクトル( σ = Hc )の表式は，
σ Iα I β = ∑ H Iα I β , Jα J β cJα J β( Jα J β )β= ∑ ∑ u pq  Iα E αpq Jα δ I β J β + I β E pq
J β δ Iα Jα  cJα J β( Jα J β ) pq+

1
( pq | rs )  Iα Eαpq Ersα Jα δ I β J β + I β E pqβ Ersβ J β δ Iα Jα  cJα J β∑∑
2 ( Jα J β ) pqrs

+ ∑ ∑ ( pq | rs ) Iα E αpq Jα( Jα J β ) pqrs

(1)β
I β E pqJ β cJα J β

1 core(2)∑ ( pr | rq ) ）2 rで与えられる．MC-SCF 法においては，必要な 2 電子積分はすべてメモリー上に収納できること，また，
coupling 係数 Iα Eαpq Jα は string Jα に対する 1 軌道置換のテーブル（ Iα Eαpq Ersα Jα（ただし， u pq = hpq −は 1 軌道置換を 2 回繰り返すことに対応する）であるから，104 次元以下程度の string 空間に対しては，これらもメモリー上に収納できることから，string をもとにすると，すべて，メモ
リー上で高速に計算できることがわかる．しかしながら，一方では，(Jα Jβ)は Jα と Jβ に対する和が独立でなく Jα と Jβ の組について和をとることを示し，それゆえ，一般の MC-SCF の場合には，Kronecker の δ があっても，CAS の場合のように和は簡約されず，二重和のまま残ってしまう．
これまで，以下のような MC-SCF 法を試行してきた．
１．Quasi-complete active space (QCAS) SCF 法
２．Parent configuration plus singles and doubles (P+SD) SCF 法
３．Occupation specified space (OSS) SCF 法
４．String product space (SPS) SCF 法
１の QCAS は CAS の積空間として定義され，size-consistency，SCF のよい収束性など，CASの多くの特長を保持している．また，２の(P+SD)空間は，少数の親配置とその 1, 2 電子励起
配置からなり，他の空間と比べて小さい次元で極めてコンパクトに波動関数を表現することができる．３の OSS は，QCAS の一般化にあたる空間で，QCAS と同様に積空間として定義

されるが，OSS を構成する空間のそれぞれは CAS である必要はない．size-consistency, 収束性の良さを失うが自由度はより増している．
４の string product space (SPS)は，昨年度の討論会でも発表したように[1]，α-string 空間 Aと β-string 空間 B との積空間として定義される空間である．

P =A⊗B(3)

一般にA とB は complete ではないので，得られた波動関数は，size-consistency を満たさず，また，各軌道配置に対してすべてのスピン関数を含むわけではないので，S2 の固有関数である保証もない．波動関数の変分空間としては一種奇妙な空間であると言える．しかし，それでもなお SPS はいくつかの魅力的な特長を有している．一つは，α-string 空間A と β-string 空間B が独立であるため(1)式の和 Σ は独立な和 ΣΣ となり，最後の項を除いて一重和 Σ ( Σ )( Jα J β )JαJα J βJβに簡約され，CAS の場合と同様の効率のよい表式に帰着することであり，また，もう一つは，多重励起を部分的に含むので，例えば，結合の解離のよい記述が可能なことである．
ここでは，三重結合の解離曲線の結 E / hartree
果について示す[2]．三重結合の解離 -112.4
を記述するには，異核の場合には 3電子励起，等核の場合には 6 電子励起 -112.5が必要である．三重結合の異核 2 原子-112.6分子の例として CO 分子を，三重結合の等核 2 原子分子の例として N2 分子 -112.7を採用した．
GMC(SD)-SCF
図１は，CO 分子について，8 電子 -112.8
SPS-SCF13 軌 道 を 用 い て 構 成 し た SPS[Hartree-Fock(HF) 配 置 を 親 配 置 と し -112.9
CAS-SCFて string 内で 1,2 電子励起], GMC(SD)
空間 [HF(P)+SD], および，CAS によ -113.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
る MC-SCF 法の結果を比べたものでr(C-O) / Å
ある．GMC(SD)-SCF が解離極限では
図１
大きく外れるのに対し，SPS-SCF は

CAS-SCF の結果をよく再現していることがわかる．
図２は，
N2 分子について，
6 電子 12 軌道を用いて構成した SPS [HF 配置を親配置として string
内で 1,2 電子励起], SPS(L) [HF, (HOMO)2→(LUMO)2, (HOMO)2→(2nd LUMO)2 を親配置], および，CAS による MC-SCF 法の結果をE / hartree
同様に比べたものである．SPS-SCF は-108.7
5, 6 電子励起を含まないため，解離極限を正しく記述しない．これが，-108.8SPS(L)-SCF では，7 つのうち 2 つの解離配置を親配置に含めることによっ-108.9てほぼ完全に回復していることがわ
SPS-SCFかる．
-109.0
その他の MC-SCF 法 1-3 の例についSPS(L)-SCFては当日発表する．
-109.1[1] 白井，中野，平尾，分子構造総合
討論会、3C13(2002)，神戸[2] S. Shirai, H. Nakano, and K. Hirao,
in preparation

CAS-SCF-109.2
0.5

1.0

1.5

2.0

r(N-N) / Å
図２

2.5

3.0

3.5

2Ta02露に相関した電子状態理論の発展（名大院情報科学）天能精一郎

量子化学では多くの研究グループが より大きな系を より精密に 計算可能な手法を確立し，多くの実験研究を理論的手法で置き換えるという目標に向かって開発競争を進めている。次世代
分子軌道理論で最も重要なキーワードはリニアスケーリング(Linear Scaling)と露に相関した(Explicitly correlated) 手法であり，今後十年間は非経験的分子軌道法で大街路となると予想されている。ここでは後者の発展を総括的に報告する。

1. 緒言
高速計算機の発達に伴い，高精度分子軌道理論の開発が 1980 年頃から精力的に行われてきた。
結合クラスター理論やその応答理論，多参照配置間相互作用法が多くのプログラムパッケージに採用され，高い信頼度で幾何構造や励起エネルギー，反応熱を計算する事が可能となっている。
一方，これらの手法で本質的である電子の短距離相関の取り扱いは，一電子基底関数の角運動量子数に対する収束性が ( L + 1)− 3 と極めて遅い事が知られている。最近の研究から，頻繁に使用されている二倍分極基底関数では二〜三原子分子の原子化エネルギーですら 10-30kcal/mol の系統的な基底関数打ち切り誤差を生じ，化学的精度と呼ばれる１kcal/mol 程度の誤差内で実験を再現する
ためには cc-pCV5Z や 6Z といった十倍規模の基底関数系を用いなければならない事が分かってきた。このような理論計算が現実に可能なのは数原子分子に限られ，応用量子化学の分野では分子軌道理論の真価が発揮しきれていないというのが現状である。
2. 露に相関した理論
2.1 カスプ条件とジェミナル展開
二体相関が本質的な動的電子相関を一電子基底関数の積和で近似する事は明らかに非効率的である。これを直接的に改善するのが露に相関した二体関数（ジェミナル）の導入である。電子相
関に対する近距離での漸近的振る舞いは，一重項と三重項対についてそれぞれ一次摂動から導かれる s−波と p−波のカスプ条件として知られている[1]。
Ψ (1) (x1 , x2 ) = [r12+ O(r122 )]Ψ (0) (x1 , x 2 )2( j + 1)(1)

ここで j は相対座標を極座標で表した時の角運動量子数であり，一重項と三重項対について 0 と 1の値を取る。この電子間距離に線形な相関は原点で滑らかでないので，占有対関数に短距離で電子間距離を乗じた二体基底関数を加えれば CI 展開の収束性を著しく改善出来る。古くから小さい
原子・分子の精密計算には Hylleraas CI[2]や露に相関したガウス型ジェミナル(GTG)[3]を用いた多体摂動論が用いられてきたが，三電子や四電子積分のために実際の化学の問題の応用される事が困難であった。変分モンテカルロ法ではこれらを数値的に取り扱う事で困難さを避けている。ジェミナルで相似変換された有効ハミルトニアンを用いる Transcorrelated 法[4]も相互作用が三体で閉じるというメリットがあるが，得られる解は変分的ではない。
2.2 多電子積分の取り扱い
大きな技術革新は R12 理論で導入された RI(the Resolution of identity)[5]である。RI では三電子や四電子積分を二電子積分の積和で展開する。
pqr X (r12) Y (r13 ) stu ≅ ∑ pq X (r12 ) vt v r Y (r12 ) su

(2)v

これにより，露に相関した理論の計算時間は通常の一電子基底関数展開と同等のスケーリングとなり，幅広い応用が可能な手法となった。中間関数 v の角運動量子数は二次の摂動論の場合，3LOCC

で収束する。ここで LOCC は占有軌道の角運動量子数である。従って，遷移金属の d-殻を含むような系の計算には i-殻までの中間関数が形式上必要である。我々は三中心積分を用いた密度フィッ
トと RI を組み合わせた以下の標識を提案した。
pqr X (r12 )Y (r13 ) stu ≅ ∑ p A s Aq X (r12 ) vt vr Y (r1 2 ) 1u

(3)

Av

密度フィットは当初 Manby によって R12 法に適用されていたが[6]，高速化の目的であり収束性は改善されていない。新たな標識では中間関数の角運動量子数が 2 LOCC で収束し，著しく精度が高められる事が示された[7]。密度フィットについては今後も発展が期待される。
2.3 スピンコンタミネーション
量子モンテカルロ法や Transcorrelated 法で用いられる相関因子 C = exp[ G ] の推進演算子はスピン空間座標に依存した形をしている。

G = ∑ g ( x i ,x j )

(4)

i> j

これらの方法で用いられてきた推進演算子 G は反平行スピンの一重項と三重項対関数を区別しな
いので二つのカスプ条件を同時に満たす事が出来ず，スピンコンタミネートした解が最安定となる。最近我々は，空間座標を交換する演算子 pij を含むスピンフリーの推進演算子

G = ∑[ g (d ) (ri ,r j ) + g ( x) (ri , r j )pij ]

(5)

i> j

が二つのカスプ条件を完全に満たす事を示した[8]。カスプ条件から各ジェミナルは漸近的に
3
1
g (d ) (r1 ,r2 ) = r12 + O (r122 ) , g (x ) (r1 ,r2 ) = r12 + O(r122 ) ,
8
8

(6)

となる。これを短距離ジェミナルと組み合わせる事により露に相関した方法の理論構成が簡素化され，結合クラスター理論や多参照 CI 法などへの適用が著しく容易になる。

Pair correlation energy / mEh

3. 数値的な有用性
一例として(21s15p12d8f4g)基底を用いた一価の銅原子に対する二次摂動対相関エネルギーを通0
常の MP2 値からの差として示す。短距離ジェミナルを用いた場合，通常の基底関数展開ではかなり大規-5
模な基 底 関 数 を 用い て も再 現 出 来 な か っ た 約-10100kcal/mol（150mEh）の相関エネルギーを取り込む-15
事が可能となっている。当日は，多参照配置間相互-20
作用法に適用した反応熱や励起エネルギーの結果も合わせて講演を行う。
-25
-30

MP2 :

-1477.60 mEh

3p2
2p3d
3p3d
3d2

2s3d
3s3d
2p2
2p3p

2s2p
2s3p
3s2p
3s3p
1s3d

2s3s
3s2
1s2p
1s3p

1s2
1s2s
1s3s
2s2

MP2(Geminal) : -1628.05 mEh[1] T. Kato, Commun. Pure Appl. Math. 10 151 (1957); R.
-35
T. Pack and W. Byers-Brown, J. Chem. Phys. 45 556-40(1966); W. Kutzelnigg and J. D. Morgan III, J. Chem
Phys. 96 4484 (1992) [2] E. A. Hyleraas, Z. Phys. 54 347(1929) [3] K. B. Wenzel, J. G. Zabolitzky, K. Szalewicz, B. Jeziorski and H. J. Monkhorst, J. Chem. Phys.
85 3965 (1986) [4] S. F. Boys and N. C. Handy, Poc. Roy. Soc. A. 310 43 (1969); S. Ten-no and O. Hino,
International Journal of Molecular Sciences 3 459 (2002), and references therein. [5] W. Kutzelnigg, W.
Klopper, J. Chem. Phys. 94 1985 (1991) [6] F. R. Manby, J. Chem. Phys. in press (2003) [7] S. Ten-no and
F. R. Manby, J. Chem. Phys. in press (2003) [8] S. Ten-no to be submitted.

2Ta03スピン軌道 CI 法による簡単な分子の非断熱過程の理論的研究（慶大理工） 藪下 聡

【序】量子化学計算で考慮する相互作用は、多くの場合、静的な電荷粒子間に働くCoulomb力のみで、磁気的な相互作用は無視している。しかし電子は単に負電荷をもった静的な粒子ではなく、その固有のスピンに基づく磁気モーメントをもつ動的な粒子である。このため現実の分子系では他の様々な相互作用が含まれ、多彩な化学現象を引き起こす。その中でも、特にスピン軌道（SO）相互作用は、分子内磁場と電子スピンの相互作用を表現するもので、例えば、一重項と三重項の励起状態間の項間交差（ポテン
シャル面の非断熱相互作用）や、光励起による一重項基底状態から三重項励起状態への禁制遷移の観測など、いわゆる「重原子効果」の原因としてよく知られている。特に重原子のSO相互作用は、局所的な非断熱相互作用の原因としてだけでなく、解離領域のPESの概形を決定するほどに重要である。１)
SO相互作用を理論化学の立場から考えると、これは単に相対論効果の一つである。我々は、その多種多様な計算手法の中で、実際の問題に応用することを念頭におき、一成分のSOCIプログラムの開発・応用を行ってきた。比較的簡単な分子であっても、電子相関とSO相互作用に起因する配置混合は、複雑である場合が多い。計算時間の短縮のために、また波動関数の解析を容易にするために、スピン依存
ユニタリー群、二重群、時間反転対称性などを駆使している。２)
SO 相互作用は、原子が重くなるにしたがって強くなる。しかしいつも重原子系の SO 相互作用が「顕著」で「面白く」て、軽原子系では「つまらない」かというと、それも極論である。確かに I 原子の SO 分裂は 7603 cm-1 = 0.94eV と非常に大きく、開殻電子状態では非常に顕著になる。
それに対し、Cl 原子の SO 分裂は 881 cm-1 = 0.11eV と、Ｉ原子の 1/10 で、無視されることが多い。しかしいくら小さくても有限な値を持つので、分子の解離過程に何らかの影響を与えるはず
で、小さいからこそ顕著に効く場合もある。実際、Band, Freed, Kouri らは、分解生成物が有限の
電子角運動量を持つ場合、解離領域でその角運動量の組み替えが起こり、必ず Born-Oppenheimer近似が破れることを明確にしている。本発表では、そのような例として、ハロゲン分子の光分解反応への応用と、また、O2 分子の Herzberg 吸収帯における回転線の微細構造の研究例を発表する。
【結果】ハロゲン分子の光分解反応
ハロゲン分子 X2(X=Cl,Br,I) の光分解反応は、古-29.75て X(2P3/2)と X*(2P1/2) がともに生成するが、その分岐機構や、J=1/2, 3/2 の角運動量の mJ 分布など、最近、様々の進展が見られる。特に興味深い点は、 同
じ SO 準位に相関する異なる状態間の非断熱遷移まで詳細が明らかになったことである。図１に Cl2 のポテ

Energy / Hartree

くから研究されてきた。その第一吸収帯の分解によっ3rd Ω=1u (3Σ1u + (1441))-29.76-29.77

Cl +ClA Π1uと記す。最近、λ = 355 nm の光分解3において、ともに X + X の解離極限に相関する CΠu (1u(2))→ A 3Π1u (1u(1))の非断熱遷移（図 1）が、1

B 3Π0+u3-29.79状態を 1

Cl +Cl*-29.78

ンシャル曲線を示す。ここで n 番目の Ω = 1u 断熱(n)u

Cl*+Cl*2nd Ω=1u (C 1Πu)Req
4

X 1Σ g+
5

6

7

8

9

10

R / Bohr

図１ Cl2 のポテンシャル曲線

Cl2 では起こるが Br2 では起こらないことが見出された（Kitsopoulos）。また、Cl2 では C 1Πu (1u(2))→ 1u(3) や（松見、川崎）、1u(3)→ 1u(4) の（Zare）非断熱遷移も観測されている。これらの非断熱遷移の原因は、概ね、その波動関数の角運動量の結合様式が、分子領域における LS 結合から解離領域における jj 結合に変化するためで、交換相互作用と SO 相互作用の拮抗である。図２に、SOCI 波動関数から求めた、非断熱
結合要素、gkn = < 1u( k ) |d/dR| 1u(n ) > のうち、g21, g34 を示す。
要するに角運動量の結合様式の変化に基づく分子波動関数から原子波動関数への変化を表している。我々は、これらの情報をもとに、半古典的に断熱電子状態の変化を調べ、その非断熱性を詳細に解析した。3)

図 2 Cl2 の gijCI の R 依存性

その結果、(i) C Πu (1u )→ A Π1u (1u )は非交差型で、その遷移確率は X が重くなるにつれて小さくなる。
(2)

1

(1)

3

(ii) 逆に、1u(3)→ 1u(4) は交差形で、その遷移確率は X が重くなるにつれて大きくなる。
(iii) また、C 1Πu (1u(2))→ 1u(3)が実質的に観測されるのは、X=F,Cl だけで、X=Br,I では無視できることなどが明らかになった。また、その主な原因は、(i) (ii) に対しては g21, g34 の SO 相互作用の X 依存性に、また (iii) に対しては、X=Br,I では遷移の起こる断熱エネルギー差（X の SO 分裂に対応する。図１参照）が大き過ぎるためである。
酸素分子の Herzberg 吸収帯の吸収確率
O2 の Herzberg 吸収帯はオゾン生成の key step であり、古くから多くの実験がなされてきた。3
つの禁制遷移からなり、その大部分は X(3Σ‑g)→A(3Σ+u)の第一吸収帯である。いずれも、双極禁制であるため、近年までその強度の正確な解析はなかった。Yoshino らによって回転の微細構造を含む吸収スペクトルが観測され、また England らにより、SO 相互作用と分子回転に基づくコリオ
リ相互作用（L‑uncoupling）を用いた intensity borrowing 機構（図 3）が考えられた。我々は、この単純な摂動論的モデルで回転線の吸収強度を計算しても必ずしも実験値との一致は良くない
が、正確な１次の摂動波動関数、あるいは、変分的な SOCI 波動関数を用いて、遷移モーメントを計算することにより、良い一致を得た（図４）。4)このように delicate な遷移モーメントの計算には、
丁寧な理論的考察が必要である。

SO

Z
2 00

X 3Σ-g3ΠgSO

X 3Σ-g

+

µ⊥

+

3Πg

L-uncouplingA 3Σ uA 3Σ uµ⊥

図３ intensity borrowing モデル

X
M

-26

2

||

+

A 3Σ u-1µB 3Σ-uIntegrated Line Strenght (10 cm cm )X 3Σ-g

1 50

o

P 12

o

P 1 2 (e x p )

o

P 23

o

P 2 3 (e x p )

Q

P 32

Q

P 32( e x p )

1 00

50

0

0

5

10

15

20

25

30

35

N ''

図.4 理論的回転線（ v' = 7 ← v' ' = 0 ）本研究は、浅野由花子君、竹上竜太君との共同研究によるものである。また、21stCOE KEIO-LCC プロジェクトに謝意。

【参考文献】1) S.Yabushita, in The Transition State - A Theoretical Approach, T. Fueno, Eds: 267, Gordon andBreach, Amsterdam (1999). 2) S. Yabushita, Z. Zhang, R. M. Pitzer., JPC A 103, 5791 (1999). 3) Y. Asano, S.
Yabushita, JPC A, 105, 9873 (2001); CPL, 372, 348 (2003); Bull. Korean Chem. Soc., 24, 703 (2003). 4) R.
Takegami, S. Yabushita, in preparation.

2Ta04金属錯体の光物性に関する電子論と電子状態シミュレーション
熊本大院自然科学

○杉本

学

１．緒言
有機エレクトロルミネッセンス素子や色素増感太陽電池などの光電変換デバイスに関する興味から、金属錯体の光物性あるいは新しい光応答機能を有する金属錯体の開発に関する研究が、近年極めて活発に行なわれている。金属錯体の光物性を理解し、制御するためには、その性質がどのような電子的因子のバランスにより実現されるのかを明らかにする必要がある。
我々のグループでは、単独で特徴的な性質を有する分子と複合化された金属錯体（ここではハ
イブリッド錯体と呼ぶ）の光物性に興味を持ち、ハイブリッド錯体の電子状態を解析する理論研究を行なっている。我々がハイブリッド錯体に注目する理由は、素性が明らかな構成分子とハイブリッド錯体の電子状態の相関に注目すれば、金属錯体の物性発現に重要な電子的因子を明確にしやすいと考えるからである。光エネルギー変換分子素子を開発する観点から、ハイブリッド錯体は活発な実験研究の対象となっており、このような相関に関する知見はハイブリッド錯体を設計する際の指導原理としても重要であると思われる。
電子的因子が複雑に絡み合った金属錯体の光物性を解析する一つの方法は、金属や配位子を置換して物性の変化を検討するものである。このような検討には、電子状態に関する計算シミュレーションは極めて有効な研究手法であると思われる。電子状態理論の進展によって電子励起状態についても定量性が向上する一方、計算機の高性能化によって現実的な分子の取り扱いが可能になっている。従って、電子状態シミュレーションを駆使して、分子構造を自由に変える “化学的”解析を計算機上で行なえば、金属錯体の光物性に関する理解が深まるものと思われる。
本講演では、我々が最近行なったハイブリッド錯体に関する電子論的研究を紹介し、電子状態計算の有用性と本計算化学研究により得られた新しい知見を報告したい。なお、以下に報告する
研究は、交換相関項を B3LYP 関数で近似した密度汎関数(DFT)法、Time Dependent DFT(TDDFT)法を用いて行なった。
２．ハイブリッド錯体で実現される電子状態の特徴
近年、“Molecular Square”等の特異な構造を有する集積型錯体が数多く合成されている。このような金属錯体は、複数の金属錯体から構成されるハイブリッド錯体と見なすことができる。
この種の錯体について、内部にある空孔が特殊な反応場となる可能性が指摘される一方、特徴的な電気化学的性質、分光学的性質を示すことが報告されている。本研究では、構成モジュールとハイブリッド錯体の電子状態の相関を明らかにする目的で、それぞれの分子の励起スペクトルと分子軌道の性質を比較検討した。
Molecular Square の例として[Pt(en)(bpy)]48+(1)[en = ethylenediamine; bpy = 4,4’-bipyridine]（図１）に関する検討を行なった。この分子の励起スペクトルを計算したところ（図 2）、4.07eV(305nm)
及び 4.35eV(285nm)に bpy π-π*遷移による強い吸収バンドが計算された。[Pt(en)(bpy)]48+ の励起状態を、[Pt(en)(bpy)]2+, [Pt(en)(bpy)2]2+, [{Pt(en)}2bpy]4+, [Pt(en)bpy]24+など、様々なモジュール錯体と比較したところ、モジュール錯体の励起状態は、[Pt(en)(bpy)]48+よりもかなり低エネルギー側にシフトすることが分かった。軌道エネルギー準位の比較を行なったところ、[Pt(en)(bpy)]48+とそのモ
ジュール錯体の電子状態の違いは、金属に配位子していない bpy の末端あるいは bpy が配位子していない Pt（II）による nonbonding orbital に起因することが分かった。同様の環状構造を有する[Pt(en)(bpy)]36+ の励起スペクトルは[Pt(en)(bpy)]48+ とよく類似しており、環状構造の形成による
nonbonding orbital の準位の変化が、[Pt(en)(bpy)]48+の電子状態を特徴付けると結論できる。

N
H2N Pt N
NH2

N Pt NH2NN
N Pt NH2
H2N(1): [Pt(en)(bpy)]4

8+

H2N
N

N
H2N Pt N
NH2

Pt

NH2
N

N
N Pt NH2
H2N(2): [Pt(en)(bpy)]3

2.0

in vacuum

305nm
305nm
f=1.609
f=1.609
π-π∗
π-π∗

285nm
285nm
f=0.928
f=0.928
π-π∗
π-π∗

305nm
305nm
f=0.073
f=0.073
MLCT
MLCT

1.0

304nm
304nm
f=0.135
f=0.135
MLCT
MLCT

311nm
311nm
f=0.087
f=0.087
d-d
d-d

280

図 1. 検討したハイブリッド錯体の構造

0.5Lowest excited state341nm(3.64eV)

0.0

6+

1.0
Oscillator Strength

H2N Pt NN

H2N

ε×10-5/M-1cm-1

3.0
NH2

300

320

340

0.0

Wavelength/nm

図 2. [Pt(en)(bpy)]48+の励起スペクトル３．ハイブリッド錯体の電子状態に対する溶媒効果
最近、Sauvage ら（J. Phys. Chem. B, 2002, 106, 6663）は、tetraphenylporphyrin(TPP)と iridium(III)bisterpyridine 錯体([Ir(terpy)2] 3+)を連結したハイブリッド錯体(TPP-Ir)を合成し、励起波長の選択によって分子内電子移動と分子内エネルギー移動を制御できることを報告した。この錯体の電子過程に関する興味から、我々はその励起状態の性質に関する検討を行なった。
孤立分子の１重項基底状態に関するポピュレーション(NBO)解析を行ったところ、TPP からIr(terpy)23+へ電子が流れこむことが分かった(Table 1)。実験的には、この分子の基底状態はアセト
ニトリル中で TPP-Ir(III)なる状態であるとされているが、我々の計算は基底状態が TPP+-Ir(II)に相当することを示唆する。この分子の 3 重項状態を計算したとこ
ろ、この状態は１重項よりも安定であり、この状態もTPP+-Ir(II)に帰属できることがわかった。

表 1. １重項基底状態における TPP-Ir の電荷分布(NBO 解析による)In vacuum
In acetonitrile

TPP
+0.88
+0.16

Ir(terpy)2
+2.12
+2.84

このような大きな電荷の偏りがある系では、溶媒和によって電子状態が大きく変化する可能性がある。このため Polarizable Continuum Model (PCM)を用いて、アセトニトリルによる溶媒効果を考慮した計算を行った。その結果、１重項状態が３重項状態よりも 1.16eV 安定となること、１重項基底状態が TPP0-Ir(III)に相当することが分かった。こ
の結果は TPP-Ir(III)なる電子状態を実現する上で、溶媒の静電的な reaction field が不可欠であることを示してい③CT(TPP→Ir(terpy)2)
TPP π→π*る。これは、この錯体で観測される電子移動とエネルギ
ー移動にも溶媒効果が本質的に重要な役割を演じることを示唆する点でも興味深い。また、我々の計算では①CT(TPP→Ir(terpy)2)②TPP(Q 帯)

1000nm 近傍の近赤外領域に吸収バンドが現れることを示す結果が得られた。実験的にはこの低い励起状態は観測されていない。今後、この吸収バンドの有無に関する実験的検証が望まれる。

図３．アセトニトリル中の TPP-Ir の励起スペクトル（計算結果）

2Ta05

機能性高分子の量子化学的設計方法と計算方法の開発(広大院理・科技団 PRESTO) 青木百合子

巨大高分子系の電子状態を効率よく計算する方法として現在までに開発してきた O(N)の計算方法である高分子の理論的合成法−Elongation 法−に、高分子の磁性・導電性・非線形光学現象などの機能設計を行う手法を組み込み、高分子の電子状態から高分子機能を予測する計算手法を開発している（図１）。さらに、電子相関の効果も含めた基底状態および励起状態の定量的な軌道相互作用解析が可能となった NBO Based CI/MP Through Space/Bond 解析法とも結合することによって、構造・物性と電子状態の関係の詳細な解析のための高分子物性・機能解析システムの構築をめざしている。今回は主に、非周期性高分子の磁性・導電性・非線形光学現象の予測に関する計算方法と現在までの計算結果について紹介する。
分子に磁性をもたせる試みが盛んに行われており、近年 S>5000 を有する有機強磁性共役ポリマーが合成された。有機化合物でありながら磁石としての性質をもつ共役系を予測するための簡単なルールを提出し、その高スピン性を CAS‑SCF 法で確認すると共に、ab initio CI/MP ThroughSpace/Bond 解析法によって、ラジカル間の相互作用の経路を調べる試みを行ってきた。これらの結果より、高スピン型有利な分子間結合のさせ方について、Non‑disjoint 性（NBMO の縮重が保たれつつ NBMO 間の混じり合いが多くなる性質）の高い超分子を与えるための結合規則−「NBMO 係
数の０の部分と係数のある＊の部分を結合させる０−＊結合が Non‑disjoint 型を与える」1,2)が有効であることを確認することができた。これを分子ユニットから高分子鎖を構築する過程に適用すると、高分子全体として強磁性としての性質を示しうる高分子構造や置換基について理論的な知見と分子設計の指針を得ることが出来る。ただし、２分子間の結合のみを厳密に調べて高分子鎖の強磁性を予測しても、全系をまるごと計算すること以外に、凝集系として実際にスピン整列が起こるのかどうかを定量性的に確認する方法が存在しない。そこで、以前より開発している
高分子の理論的重合法（Elongation 法）を非制限 Hartree‑Fock（UHF）法のレベルで計算可能とし、3〜5 重項程度の低スピン状態計算を繰り返すだけで S=>∞の状態まで効率よく計算できるよう発展させた。本方法によると、高分子の各ユニットに局在化した局在化分子軌道(LMO)が得られるが、高分子鎖を伸長する際に、取り扱う相互作用領域で、常に高スピン状態になるように結合させる。その際の Non‑disjoint 性の大きさも、得られた LMO から算出することができる。
同様の概念を分子間に適用することにより、高スピン型としての性質を発現しやすい分子クラスターおよびその分子間の配向様式に関して、理論的な予測を提示することが可能となる。結晶を構成する分子としては、不対電子スピンを有する有機分子や、ドナー分子とアクセプター分子が交互に積層した電荷移動錯体などが考えられる。しかし、結晶全体として強磁性になるか反強磁性になるかは、スタッキング方向における積層構造の状態や鎖間相互作用に大きく依存する。
分子に対して、なるべくスピンを整列させるための条件である１．非結合性軌道（NBMO）縮重系を構築する、２．NBMO 間に弱い相互作用をもたせる、を同時に満たすような積層構造の持たせ方の他、適度な鎖間相互作用をもたせることにより、系がパイエルス不安定性を持たないように制御することが重要となると思われる。擬一次元鎖を考えた場合、分子に非結合性軌道（NBMO）が存在してもそれが一次元的なカラムをなした場合、往々にして系のもつ Peierls Instability のために、２分子単位で閉殻構造になり、NBMO 間の縮重性が失われてしまう。つまり、構成分子のラジカルが凝集したときに伝導電子としての役割も加わるため、金属―絶縁化転移と強磁性としての性質は密接に関係している。これに、鎖間相互作用を与えると、Peierls Instability が抑えられ、等間隔の配列が有利になり得る 3‑8）。言い換えると、NBMO 間の縮重性が保たれるように分子間を等間隔にするには、２量化が抑えられる程度の鎖間相互作用が必要ということになる。
つまり、０−＊結合の考えを分子間に拡張すると、一次元鎖内および一次元鎖間において、なるべく０の部分と＊の部分が重なるようにしてスタックさせると高スピン性には有利に働く。結局、
Non‑disjoint 型でスタッキングさせ、しかも Peierls Instability が抑えられるよう適度な鎖間相互作用をもたせることが、系全体としてスピンを整列させるのに有利な条件であるとともに、導電性の発現にもつながると期待できる。これらの関係について、いくつかの系に対して本方法を適用し、高スピン性と Non‑disjoint 性の関係を検討する。
一方、Elongation 法では種々の分子を順次付加できるため、非周期的な高分子の電子状態計算に適している。構造の非周期性が系全体の物性に大きく影響する物理現象のひとつに非線形光学現象（超分極率）がある。周期性高分子であれば、周期境界条件のもとで分子（超）分極率を求める方法は開発されているが 9,10)、一般的に超分極率を発現する典型的な系として有限鎖である
Push‑Pull 有機共役高分子系がある。これは、適当な長さの一次元オリゴマーの両端に電子吸引性および電子供与性置換基をもつもので、系全体としての反転対称性の崩れに伴って超分極率を発現する。このような、分子と周期性高分子の間の中間的な大きさの領域にみられるメゾ系の電子状態を効率よく計算するためには Elongation 法が有用である。計算方法としては、最初からハ
ミルトニアンに電場の効果を導入して固有値問題を解き、全エネルギーの電場勾配を求める
Finite Field(FF) 法 と 、 電 場 の 効 果 を 摂 動 法 に よ っ て 評 価 す る Coupled Perturbed
Hartree‑Fock(CPHF)法がある。両方法を Elongation 法に組み込み中であるが、まず試験的にFinite Field 法を Mopac‑Elognation 法に組み込んだため 11)、これによる計算結果を紹介する。
本方法では、電場の効果（‑r・E ）を最初からハミルトニアンに加えて出発クラスター(Elongation
法の出発となる小さなクラスター)を計算しておき、小分子が付加するたびに、電場の摂動ハミル
トニアンを含む相互作用 Fock 行列を計算し、付加分子に近い領域に局在化した ActiveLMO を基底とした固有値問題を解く。分子を一個付加する度に、５種類の電場の大きさのもとで電子状態を計算し、数値微分によって分子分極率(α)超分極率(β,γ)を計算する。特に末端置換基効果を調べるにあたっては、置換基以外の部分まで高分子鎖の電子状態を伸長しておき、末端置換基のみを交換することによって、効率的に置換基効果を調べることができる。これら非周期的高分子の電子状態と機能に関する計算手法と現在までの計算結果について発表する。
参考文献）1)Y.Aoki and A. Imamura. 特許：特開平１１−６８２５ 「有機高分子化合物の強磁性を予測する方法」. 2) Y.Aoki and A. Imamura.
Int. J. Quantum Chem., 74, 491‑502(1999). 3)Y.
Aoki and A. Imamura, J. Chem. Phys., 103,
9726‑9737 (1995). 4)Y. Aoki, T. Tada, and A.
Imamura, Int. J. Quantum Chem., 64, 325‑336(1997). 5) T. Tada, Y. Aoki, and A. Imamura,
Synthetic Metals, 95,, 169‑177(1998). 6)T.
Tada and Y. Aoki, Int. J. Quantum Chem., 86,
401‑415(2002). 7)T. Tada and Y. Aoki, Phys.
Rev. B, 113113 1‑4(2002). 8) Y. Aoki, T. Tada
and Y. Orimoto, Phys. Rev. B, 66,193104 1‑4(2002). 9)F. L. Gu, Y. Aoki and D. M. Bishop,
J. Chem. Phys., 117, 385‑395(2002). 10)B.
Champagne, D. Jacquemin, F. L. Gu, Y. Aoki,
B. Kirtman, and D. M. Bishop, Chem. Phys.
Lett., 373, 539‑549 (2003).11) F. L. Gu, Y.
Aoki, A. Imamura, D. M. Bishop, and B. Kirtman,
Mol. Phys., 101, 1487‑1494(2003).

Elongation法
逐次的

並列的

高分子の構造決定

有限軌道基底[O(N)]

半経験的MO法

非経験的MO法

NBO based ab initio
CI/MP2 TS/TB解析

導電性高分子設計
強磁性高分子設計
高分子の非線形光学現象（超分極率の計算）（生体）高分子の協同現象の解明

機能性高分子の設計

機能性高分子の解析

2Ta06ヘム錯体の基底・励起状態とＮＭＲ化学シフト（都立大院理）波田雅彦【序】
ヘム錯体は Fe(II)と Fe(III)及び高スピンから低スピン状態までの種々の電子状態の変化を伴って生体反応に関与している。NMR スペクトルはこれらの電子状態や分子構造の変化を推測する有用な手段として用いられている。昨年、Fe(III)‑bis(cyanide)錯体において、Feに配位する CN の 13C 常磁性化学シフトが約‑2500ppm(TMS 基準)の高磁場領域に存在し、
cyanide‑imidazol(Im)錯体ではさらに高磁場シフトすることが報告された[1]。これら一連のスペクトルは heme 周辺の僅かな化学的変化を鋭敏に反映するため、heme における生体反応
メカニズムを解析する手段として注目される。
本報告は、低スピン状態の Fe(III)錯体をモデルとして、配位子やポルフィリン環の歪による電子状態の変化とC‑13 NMR スペクトルとの関連を理論的に明らかにすることを目的とし、種々の Fe(III)ポルフィリン錯体の基底・励起電子状態及び 13C 化学シフトの計算結果を報告する。本要旨では右図の L = CN‑, Im、および R = H, Ethyl(Et)基の結果を述べる。

Fig.1 cyanide-porphyrinato Fe(III)【計算方法】
基底・励起電子状態及びその電子スピン密度の計算はGaussian03 に組み込まれた SAC/SAC‑CI を用いた。スピン密度を計算するため配置選択の閾
値は十分小さな値を用いた。基底関数は Fe に(62221/6211/311) + diffuse‑sp + f、CN 配位
子には(721/41) + diffuse‑sp + d、ポルフィリン環の原子には(651/51) + d を用いた。
原子 A の化学シフト δ A は磁気遮蔽定数 σ A の差であり、 δ A = σ 0 − σ A となる。 σ 0 は TMS の
磁気遮蔽定数とした。電子スピンに依存しない σ A は Ramsey の式で求めた。この項に、常磁性分子の電子スピンに由来する Fermi‑Contact(FC)項

σ A (FC) = 4π g 2 β 2 Ψ Sˆδ A Ψ / k BTを加えた。上式は二重項の場合に成立する。本報告では、Spin‑Dipolar 項の寄与は小さいものとして無視した。
【結果と議論】(1) L = CN‑ と L = Im
表１に R = H の計算結果を示す。基底状態はいずれの錯体も(dxy)1(dxz,yz)4 が主配置であり、実験的示唆と一致する。但し、両錯体ともに(dxy)2(dxz,yz)3 配置がごく近傍(0.2 eV 以内)に存在しており、環構造や周辺の変化によって基底状態が入れ替わる可能性を示唆している。σ∗(CN)が open‑shell である配置は溶媒効果によって若干不安定化すると予想される。
化学シフトの計算値も表１に示してある。シフト値の大部分が FC 項の寄与である。

bis‑cyanide 錯体(L=CN‑)の化学シフト
では、実験値‑2516ppm の高磁場シフトを計算値‑2789ppm が良く再現している。逆に、2A2u の化学シフトは異常に低磁場シフトしている。このようなシフト傾向の大Table 1. Energies and Chemical Shifts of bis‑
cyanide and cyanide‑imidazol complexes.
Shift
(ppm)
bis(cyanide) porphyrinate Fe (L=CN‑, R=H)

Term Energy(eV)Main Configuration

2Eg

0.00

(dxy)2(dxz,yz)3

‑2789

2A2u

0.09

d6[σ(CN)‑π(Por)]1

16557

0.57

(dxy)1(dxz,yz)4

‑3502

2A1u

0.74

d6[s(CN)‑π(Por)]1

‑242

1A1g

1.2

Fe(II) (closed‑shell)

171

きな相違は、下図に示すように、CN 配位
2B2g
N
N
A2u
B2u N
Eg
C

C

C

Fe

Fe

Fe

C
N

C
N

C
N

cyanide‑imizazol porphyrinato Fe (L=Im, R=H)子の炭素上の電子スピン密度がスピン分
2

A

0.00

d6[σ*(CN)‑4p(Fe)]1

10870

2

A

0.13

(dxy)2(dxz,yz)3

‑3312

2
体(L=Im)では、(dxy)1(dxz,yz)4 と(dxy)2 A

0.19(dxy)2(dxz,yz)3

‑3003

A

0.46

(dxy)1(dxz,yz)4

‑4026

2
る 。 し か し 、 基 底 状 態 で な い (dxy)1 A

0.49

d6[p(Por)]1

極かスピン非局在化のどちらで生成したかに依存している。cyanide‑imidazol 錯(dxz,yz)3 のいずれも観測の傾向を再現す

2(dxz,yz)4 配置の化学シフトが観測値と最も良く一致する。この点はさらに理論・実験両者の検討が必要である。
(2) L=CN‑, R=Et
ポルフィリン環の meso 位水素を Ethyl 基で置換すると環が歪み、D4h 対称性は湾曲した S4となる。環の湾曲により dxy 軌道は不安定化
し、実験で示唆されているように‑932

Table 2. Energies and Chemical Shifts(L=CN, R=Et).
Term Energy(eV) Main Configuration

Shift
(ppm)

2B2

0.00

(dxy)1(dxz,yz)4

4039

2E

0.04

(dxy)2(dxz,yz)3

‑3121

2B2

0.37

d6[s(CN)‑p(Por)]1

10061

2B1

1.0

d6[p(Por)]199(dxy)1(dxz,yz)4 配 置 が 基 底 状 態 と な る ( 表 ２ ) 。 但 し 、 σ∗(CN) が - 4000

Eg → E

open‑shell である配置は溶媒効果によって不安定化が予想されるの - 2000
で基底状態の議論からは除外する。
一方、化学シフトは環が僅かに歪むことによって約 7500ppm(‑3502ppm

0→+4039ppm)という異常な変化を示した。図１にポルフィリン環の歪み 2000
に対するシフト値の変化を示す。B2g と A2u の変化が相補的である理由は、 4000
環の歪みによって両電子状態のスピン密度が混合するためである。
【結論】
励起状態と常磁性化学シフトの精密な理論計算は、NMR スペクトルの

10000

14000

段となる。本研究は Fe(III)ポルフィリン錯体におけ

16000[1] H. Fujii, JACS, 124, 5936, 2002

A2u → B212000

観測結果を補い、heme 蛋白の電子状態や構造変化を解析する有力な手
る 13C 常磁性化学シフトに重要な解釈を与えた。

B2g → B2

D2h → S4
Fig. 1. 13C chemical shifts ofplanar (D2h) and ruffled (S4)complexes.

2Tp01

局所内挿法によるポテンシャル面の効率的生成（分子研、Northwestern 大化学）
○石田俊正,George C.Schatz［序］計算機性能が向上するのにともなって、電子状態を解きながら原子核を運動させていく abinitio 動力学が最近行われるようになって来た。しかし、反応動力学を正確に記述しようとするならば精度の高い ab initio 計算結果が必要であり、古典動力学で満足するとしても統計的に有意なトラジェク
トリ数を確保するのは依然として困難である。一方、大局的な解析関数によるポテンシャル面のフィットが以前から行われてきたが、複雑な関数型や豊富な経験が必要とされ、簡便とは言い難い。
そこでわれわれは局所内挿法に着目し、ab initio 法において導関数を求めなくてすむ IMLS/Shepard法を提案し、応用を行っている[1-3]。本方法では、Collins らによる修正 Shepard 内挿法[4]と IMLS(内挿移動最小二乗)法とを組み合わせている。なお、Shepard 法は 0 次の IMLS 法に対応している。ab initio
計算において導関数を求めなくてよいので、精度の高い最先端の方法論と組み合わせて用いることが可能である。さらに、反応動力学に関与しない高エネルギー領域のサンプリングを排除できること、
原子核の置換による対称性が簡単に考慮できること、精度が足りないと考えられる場合は、内挿点の追加によるポテンシャル面の精度の向上が容易なこと、当てはめに要する関数のパラメータ数が実質的に一つであることなどが利点である。
Bettens と Collins は Bayesian の定理に基づいて重み関数を決める方法を提案している[5]。比較的局
所的な範囲と比較的広い範囲を表現する２つのパラメータ(p と q)を用いた重み関数と Bayesian 解析法とを用いて、彼らはポテンシャルの精度が改善されることを報告している。そこで、この方法を
IMLS/Shepard 法と組み合わせたときの効果についても O(1D)+H2, OH+H2 系について検討した。本講演
では、Bayesian 解析の適用も含めて IMLS/Shepard 法の方法と応用について述べたい。
[方法]IMLS 法では、n 個の線形独立な基底を用いて、任意の点での値を基底の線形結合で表す。基底としては、核間距離の逆数 Z=1/R の 2 次までの多項式をとった。座標に依存する展開係数を重みつき最小二乗法により決定することにより、内挿関数が得られる。重み関数としては、

w = 1 /(|| Z − Zi ||2 +a 2 ) p で表した。具体的には、p として正の整数, a=0.03 を用いた。実際には、IMLS
法を内挿点のみに適用し、内挿点における（近似的な）二次微分までを見積もり、二次微分までを使う Collins らの修正 Shepard 内挿法を適用する。以上の手続きからわかるように、本方法では、基底展開の次数、重み関数の式が決まれば、関数のあてはめが完了し、簡便である。
[結果と考察](1)H2+H この系に対する Liu-Siegbahn-Truhlar-Horowitz(LSTH)による解析ポテンシャル
面[6]を用いて精度を検討した。ポテンシャル面から内挿点をサンプルし、その内挿によるポテンシャ
ル面ともとの解析ポテンシャル面を比較した。その際に、反応断面積などをトラジェクトリ計算により求めた。ポテンシャルのサンプリングには乱数を用いた。断面積は、LSTH 解析ポテンシャルによる断面積を 150 個程度の点の数でよく再現することがわかった。LSTH ポテンシャル自身が 267 点の ab
initio 計算結果に基づいていることを考えると IMLS/Shepard 法が優れていることがわかる。
(2)O(1D)+H2

この系は、反応の途中で H2O 分子が生成する深い井戸を持つ系で特異である。2 パラ

メータを用いた重み関数を用いると、IMLS/Shepard 法では、トラジェクトリ計算においてエネルギー
保存が悪かったので、重み関数パラメータ p=3 として計算を行った。その際、参照ポテンシャルとして Ho らによるポテンシャル面[7]を用いた。図２に反応断面積の変化を示した。点線は参照した（LSTH
解析ポテンシャルによる断面積）±（標準偏差）の誤差の範囲を示している。IMLS/Shepard 法においては、Bayesian 解析を用いても用いなくても反応断面積の±（標準偏差）への収束に要する点の数は変わらなかった。一方、修正 Shepard 法においては Bayesian 解析の効果は大きかった。示していないが、ポテンシャルエネルギーに対する rms 誤差について調べてみると、rms 誤差の改善のほうが大きかった。Bayesian 解析を用いると、rms 誤差が小さくなるが、ポテンシャル面の粗さは最終的にあまり改善されないと考えられる。
Bayesian 解析に Shepard 法のみを組み合わせた場合の rms 誤差の最小値が、IMLS/Shepard 法での最小誤差に近いことから、最良の結果どうしを比べると、Shepard 法と IMLS/Shepard 法は精度があまり変わらないことになる。しかるに、IMLS/Shepard 法では、ポテンシャル面の微分の情報を必要としないので、同精度のポテンシャル面を得るのに必要な計算量ははるかに少なくてすむ。このことは次の４原子系でも同様である。
(3) OH+H2 この系に対する Ochoa de Aspuru と Clary による解析ポテンシャル[8]を用いて、
IMLS/Shepard 法の誤差を検討した。図に、内挿点を増やしていったときのエネルギーの誤差（メジアン値からの偏差の絶対値）の変化を示す。解析ポテンシャル面から内挿点をサンプルする際には、乱数を用いた。IMLS/Shepard 法の結果を見ると、この系では、p の値にかかわらず、Bayesian 解析を用いた方が精度が向上していることがわかる。これは、O+H2 系での結果と異なっている。この系は、単純な遷移状態を持つ系に対して、O+H2 系は深い井戸を持つ系であり、そのポテンシャル面の違いがBayesian 解析を用いた結果の違いに反映していると考えられる。次に、重み関数を制御するパラメー
タ p について見ると、p=3,6,9 のうちでは、p=6 の結果が最も誤差が小さいが、この計算では、内挿点の数を数千取ったためであると考えられる。内挿点が多くなると、ある内挿点に着目した場合の周囲の内挿点が増え、近くの点のみの値を使う方が精度が向上することを反映している。Shepard 法と比較すると、この系の場合、内挿点数が同じ場合では Shepard 法の精度がよい。ただし、Shepard 法では、
内挿点における２次微分までの情報を用いているので、数値計算でそれを賄おうとすると、6×7/2=21
倍程度の点が必要になる。そのことを考慮すると、IMLS/Shepard 法の方が優れていると考えられる。
講演では、５次元以上への拡張についても述べたい。

図１ O(1D)+H2 系における断面積変化

図 2 OH+H2 系におけるエネルギーの誤差変化[参考文献] [1] T. Ishida and G. C. Schatz, CPL 314, 369(1999). [2] T. Ishida and G. C. Schatz, J. Comput. Chem. 24(9), 1077-1086 (2003).
[3] K. H. Kim, Y. S. Lee, G.-H. Jeung, and T. Ishida, J. Chem. Phys. in press. [4] T. Ischtwan and M. A. Collins, JCP 100, 8080(1994). [5] R.
P. A. Bettens and M. A. Collins, JCP 111 816 (1999) [6] D. G. Truhlar and Horowitz, JCP 68 (1978) 2466. [7] Ho et al. J Chem Phys 1996,
105, 10472. [8] G. Ochoa de Aspuru and and D. C. Clary, J. Phys. Chem. A102, 9631(1998).

2Tp02

自由エネルギー
自由エネルギー勾配法
エネルギー勾配法を
勾配法を用いた凝縮相
いた凝縮相ダイナミクス
凝縮相ダイナミクスのダイナミクスの研究（名大院情報科学）長岡正隆【序論】
序論】 我々が日常的に目にする化学現象の多くは、溶液中や界面で起こっている。その際、溶媒や界面は着目する反応中心に対して反応場を与える。孤立分子の安定構造や遷移状態構造が反応予測に重要な情報を提供するのと同様に、凝縮相分子の平衡構造の知見は凝縮相反応にとって重要な情報をもたらす。我々はこれまでに凝縮相分子の構造最適化法として、自由エネルギー勾配法を提案し、水溶液中分子の安定状態や遷移状態の構造最適化を実行してきた[1‑3]。本講演では、自由エネルギー勾配法の概要を説明し、その適用例から得られたいくつかの結果を紹介する。
【理論と理論と方法】
方法】

溶液反応の分子動力学(MD)シミュレーションでは、溶質分子の各原子に溶媒分子群から及ぼされる力を時々刻々計算している。この力の時間平均をとると、自由エネルギー面

G(qs)上の力
FE

S

S

S

S

F (q ) = − ∂G (q ) ∂q = − ∂VRS (q ) ∂q

S

(1)

time

s

s

が溶質分子構造 q の関数として求めることができる[2]。ここで、VRS (q ) は溶質分子内ポテンシ
ャルエネルギーと溶質溶媒間相互作用エネルギーとの和を表わす。また " time は時間平均を表わ
すが、平衡系ではアンサンブル平均" B =∫ dq (") exp(− βV ) ∫ dq exp(− βV )
B

B

(2)

B

S

に等しい。ここで、q は溶媒分子群全体の座標を表わす。また V は、溶質分子構造 q を固定した条件下で得られる溶液全体のポテンシャルエネルギーで、相互作用の評価にハイブリッド
QM/MM 法を採用すると、V = VRS (q , q ) + VMM (q ) 、 VRS (q ) = Ψ Hˆ QM + Hˆ QM/MM Ψ = VR (q ) + Ψ Hˆ QM/MM Ψ
S

B

S

B

S

(3)

と表現できる。 VMM は溶媒分子群のポテンシャルエネルギーを表わす。
sss

平衡ＭＤ計算を繰り返して溶質分子構造を更新( qi +1 = q i + ∆q i )して、或る nopt 番目の更新でゼs

ロ勾配条件 ∂VRS (q ) ∂qsk

≈ 0 が、十分精度良く満たされれば溶質分子の構造最適化が達成され安定平衡状態が求まる。このとき、力釣り合い条件∂ Ψ Hˆ QM Ψ
∂VR (qs )
=
∂qs
∂q s

=−

∂ Ψ Hˆ QM/MM Ψ

(4)

∂qs

が成り立ち、溶質自身のポテンシャル勾配と、溶媒が溶質原子に及ぼす力とが釣り合う。
ss

一方、溶質分子の構造 q i から q i+1 への最適化過程に伴う自由エネルギー変化 ∆Gi は、溶質構s

造 qsi の下でアンサンブル平均を用いて統計摂動法[4,5]により与えられる。こうして q 0 から最適化
構造 q n までの安定化自由エネルギーは、 ∆Gstabi = ∑ i opt=0 ∆Gi で得られる。
n

S

opt

i

【結果と
結果と考察】考察】A. 水溶液中アンモニア
水溶液中アンモニア‑
アンモニア‑水分子対の
水分子対の構造最適化 [6]水溶液中のアンモニア‑水分子対（H3N…H2O）の構造最適化を実行した。QM/MM法としては、Ruiz-López等の方法１[7]に従って、H3N…H2O対にはPM3ハミルトニアンを、溶媒水分子にはTIP3Pモデルを用いた。図１に、最適化ステップ毎の二乗平均変位を示した。第８ステップまでに分子対構造が最適化されている様子がわかる。図２には最適化構造を示した。N原子と水素結合していないH原子を含む

OH結合が孤立分子対よりも長くなり、溶媒水系のミクロ構造を考慮した結果生じた現象と解釈できる。これは簡単な誘電体モデルでは再現できない。
-2

RMS displacement [Å]2.0x10

1.5

0.998

1.0

107.2

0.980

1.812

0.5

1.001

178.9

107.4

0.964

0.0
0

2
4
6
Optimization Step Number8

図１．最適化ステップ毎の二乗平均変位

図２．H3N…H2O対の安定状態構造（結合長:Å、結合角:度）

B. NH3 + CH3Cl → CH3NH3+ + Cl– (遷移状態完全最適化の遷移状態完全最適化の最初の最初の例) [8]自由エネルギー勾配法を用いて、水溶液中メンシュトキン反応 NH3+CH3Cl→H3NCH3++Cl-に対する遷移状態構造を求めることに成功した(図 3)[8]。また気相中に比べて水溶液中の反応経路では遷移状態が反応系よりに大きくずれ、より電荷分離するように構造変形することが示された(図 4)。

(a)
1.10
1.97

91.2
2.07

1.00
109.4

(b)
1.11
1.66
1.01

80.2
2.24

110.4

図３．メンシュトキン反応の(a)安定状態構造と(b)遷移状態構造（結合長:Å、結合角:度）

図４．気相と水溶液中におけるメンシュト+
キン反応NH3+CH3Cl→H3NCH3 +Clの自由エネルギープロファイル【参考文献】
参考文献】[1] N. Okuyama-Yoshida, M. Nagaoka and T. Yamabe, Int. J. Quantum Chem. 70 (1998) 95.
[2] M. Nagaoka, N. Okuyama-Yoshida and T. Yamabe, J. Phys. Chem. A102 (1998) 8202.
[3] N. Okuyama-Yoshida, K. Kataoka, M. Nagaoka and T. Yamabe, J. Chem. Phys. 113 (2000) 3519.
[4] R. W. Zwanzig, J. Chem. Phys. 22 (1954) 1420.
[5] U. C. Singh, F. K. Brown, P. A. Bash and P. A. Kollman, J. Am. Chem. Soc. 109 (1987) 1607.
[6] Y. Nagae, Y. Oishi, N. Naruse and M. Nagaoka, J. Chem. Phys, 印刷中.
[7] F. J. Luque, N. Reuter, A. Cartier and M. F. Ruiz-López, J. Phys. Chem. A104 (2000) 10923.
[8] H. Hirao, Y. Nagae and M. Nagaoka, Chem. Phys. Lett. 348 (2001) 350.

2Tp03電子・核波束動力学法の開発と強光子場化学への応用(東北大院理) 河野裕彦

【序】高強度レーザーは，原子・分子科学や高エネルギー宇宙物理などの基礎科学に新たな展開をもたらすのみならず，エネルギー変換や物質加工への応用面も有し，きわめて学際的・総合的な科学を構築する手段として認識されつつある

1)．その出現は光と分子の相互作用の研究を非共鳴領域へと拡大させた．電子的に非共鳴な低振動数(~800nm)の強い光は無極性分子を偏光方向に配列させるなど新たなテクノロジーの手段としても注目を集めている．クー
ロン力と拮抗する強度領域（ I > 1013 W/cm2 ）に達すると（水素原子の 1s 軌道の電子が原子核から受ける引力は光強度 3.5 × 10

162

W/cm に相当），分子の電子波動関数は光の半周期（~1fs）以内でも劇的に変化し，トンネル型イオン化など数々の非摂動論的現象が起きる．
クーロン爆発時の分子構造の同定が可能になってきた現在

2)，紫外域高強度光源による極短時間イオン化による化学反応の実空間追跡も視野に入ってきた．分子は強レーザー場中で様々な様相のダイナミクスを呈し，それ故，このような強光子場化学とも呼ぶべき研究領域
3)

は，いわゆるレーザーを使った「化学反応の量子制御」 と表裏一体の関係にある．我々は，「分子の中をレーザー電場 ε (t ) によって電子がどのように動くか」を知ることが不可欠であると考え， ε (t ) と相互作用する分子の電子と核の波束（時間変化する波動関数）ダイナミクス
を解く計算法の開発に携わってきた．現在， H2 のような少数粒子系に対しては、その電子と核の全系に対する時間依存 Schrödinger 方程式を厳密に解くことが可能になってきた．ここでは、電子・核波束動力学の立場から得られた強光子場化学の研究成果を概説し、最先端レーザー技術がもたらす新たな化学の展開における理論化学の役割についても議論する．

【電子波束計算法の開発と強レーザー場中の分子の電子ダイナミクス】 クーロン系の波束を
4)

グリッド空間表示で精度よく時間発展させる Dual Transformation 法を開発した ．まず，フェムト秒の強レーザー場と相互作用する H 2+ に適用し，その電子も核も含めた全クーロン系の波5)束計算を Born-Oppenheimer 近似など一切の近似無しに行った（ab initio 動力学） ．さらに，２電子系分子の電子波束ダイナミクスに適用し，強レーザー場中の２電子の相関を時間領域
6)で見ることに成功した ． H2 では，光サイクルの半周期の間に電子が核間を移動し，局在イ
オン結合性状態 H+ H− と H− H+ が交互に生成し，不安定なイオン結合性状態からトンネルイオ
ン化することを明らかにした．つまり，強レーザー場中では実効電荷の大きい電子対が容易に生成し，局所励起よりも大域的な核間の電子移動が起こる．このような電子の動きはレー
ザー電場 ε (t ) の時間変化に追従する時間依存断熱電子状態（ ε (t ) との相互作用も含めた瞬間的な電子ハミルトニアンの固有関数）とそれらの間の非断熱遷移の考えによって記述できる．
全系の波動関数を時間依存断熱状態に射影することによって，トンネルイオン化が平衡核間距離よりずっと長い核間距離で促進される現象（Enhanced Ionization と呼ばれ，クーロン爆発を誘発する）を定量的に解明した

6,7)

．

現在では，自然軌道の時間発展に基づいた新しい Multi-reference time-dependent Hartree-Fock法を開発し，多原子分子の多電子ダイナミクスを計算することができつつある．この方法により，多原子分子のイオン化過程や反応過程を支配している強レーザー電界誘起分子内電子移動のダイナミクスをより詳細に記述することができるであろう．
【強電場中の分子の電子状態計算と強光子場化学】 H2 などの厳密波束計算から，強い光の場

の中でも核が時間依存断熱電子状態のポテンシャル面上を動くという考えが成立することを
示した

5-7)

．分子軌道法に電場を取り込むことによって時間依存断熱状態を求めることができ

8)，さらに，時間依存断熱状態間のレーザー電界誘起非断熱遷移確率も計算することができる

9)．このような手法で，強光子場中の多原子分子のダイナミクスの特徴（多価イオン生成や分子構造変化）を明らかにした．例えば,電子的に非共鳴な 800nm の光と相互作用する CO 2 では，２価カチオンの段階で，屈曲しながら２つの C-O 結合が同時対称的に伸びる過程が存在9)することを見出した ．強度を上げると，一つの C-O 結合だけが切れる確率は減少し，強レ
ーザー場中での反応制御の可能性を理論的に明らかにした．このような特異な反応ダイナミ
クスは，実験的に同定されたクーロン爆発する瞬間の CO 2 多価カチオンの変形構造を説明す
る

10)

．

また，エタノールの選択的解離反応も解析した．中性では C-C 結合と C-O 結合はほぼ等し
い結合エネルギーを持つが，山内・神成ら

11)，C-C 結合の方が C-O 結合より切れやすいという実験結果を報告している．計算結果によれば，中性のエタノールが構造変化する前に，まずイオン化が起こる．生成した１価カチオンの段階で，C-C 結合が選択的に解離するという結果を得た．最低断熱状態とより高い断熱状態の時間依存ポテンシャルが短い C-C 結合距離で交差し，２つの状態間の非断熱遷移により解離が促進されている．
12)【制御】Levis らは高強度の波形整形パルスによってアセトフェノンの組み替え反応の収率を大きく変えることに実験的に成功した．強レーザー場で誘起された電子運動によって核が動き出す．このような電子−核相関を光で“調整”することによって反応を制御することができる．今後は，時間依存断熱電子状態法をそのような強レーザー場による化学反応制御の問題にも適用していきたい．また，我々は，電子的に共鳴した波長領域の位相制御された２パ
ルスの時間差を変えると，分子内電子移動ダイナミクスやイオン化を精密に制御できることを見出している．アト秒時間スケール（ 10−18 s）の電子運動の理解に基づいた化学反応の極
短時間・局所空間制御ももはや夢ではない．
「電子運動の法則」の確立は，電子状態の理論に動力学の視点を導入するという意味に限っても理論化学に新たな潮流をもたらすものと考えられる．また， 「電子・核量子波束動力学」は，現在開発が進行しているアト秒パルス

13)が拓くフロンティア「分子中の電子運動の実時間追跡とその化学反応への応用」の中心的な理論として，その重要性はさらに増していくはずである．
1) K. Yamanouchi, Science 295, 1659 (2002);

2) H. Hasegawa et al., Chem. Phys. Lett. 349, 57(2001);

3) Y. Ohtsuki, M. Sugawara, H. Kono, and Y. Fujimura et al., Bull. Chem. Soc. Jpn. 74,1167 (2001);

4) H.

Kono et al., J. Comput. Phys. 130, 148 (1997); I. Kawata and H. Kono, J. Chem. Phys. 111, 9498 (1999);
5) I. Kawata, H. Kono, and Y. Fujimura, J. Chem. Phys. 110,11152(1999); 6) K. Harumiya, I. Kawata, H.
Kono, and Y. Fujimura, J. Chem. Phys. 113, 8953(2000); Phys. Rev. A 66, 043403 (2002); 7) H. Kono, Y.
Sato, Y. Fujimura, and I. Kawata, Laser Phys. 13, 883 (2003);(2001);

8) Kono et al. J. Phys. Chem. 105, 5627

9) Y. Sato, H. Kono, S. Koseki, and Y. Fujimura, J. Am. Chem. Soc. 125, 8019 (2003); 10) A.

Hishikawa, A. Iwamae, and K. Yamanouchi, Phys. Rev. Lett. 83, 1127 (1999); 11) R. Itakura et al., J.
Chem. Phys. in press.;511 (2001).

12) R. J. Levis et al., Science 292, 712 (2001). 13) M. Hentschel et al., Nature 416,

2Tp04SAC‑
SAC‑CI 法による精密理論スペクトロスコピー（京大院工）江原 正博【序】SAC‑CI 法は分子の様々な電子状態を精度良く記述する理論として確立しており、多くの化学現象に応用されてきた。最近、SAC‑CI プログラムは Gaussian03 を通じて公開され、今後、励起状態の化学の研究に多いに活用されることが期待される。本講演では、SAC‑CI 法の応用が期待される一つの分野である「精密理論スペクトロスコピー」に関する最近の研究から、多電子過程の状態の研究、解析的エネルギー微分法による励起状態の構造やダイナミクスの研究、内殻励起・内殻イオン化スペクトルの定量的な帰属や構造緩和に関する研究を報告する。
(1) 価電子イオン化スペクトルのサテライトピーク
価電子イオン化スペクトルのサテライトピーク
SAC‑CI general‑R 法[1]は、多電子過程を精密に記述する理論である。価電子イオン化スペクトルに観測されるサ
テライトピークは多くの shake‑up 状態により構成され、それらの状態は多電子過程で記述される。これまでにgeneral‑R 法を多くの系に応用し[2]、スペクトルの詳細を明らかにしてきた。一例として、図１に N2O のイオン化ス
ペクトル[3]を示す。実験スペクトルでは、相関ピークによる５つのバンドが観測されているが、これらのスペクトルの形状を理論は詳細に再現している。ピーク I,Ⅱはそれぞれ 2Πおよび 2Σ状態のサテライトピークであり、ピークⅢ〜
Ⅴは(5σ)‑1 および(4σ)‑1 状態が分裂した相関ピークである
図 1. N2O 分子の価電子イオン化スペクト
ル(a) dipole (e,2e)と(e,2e) と(b) SAC‑
SAC ‑ CI

ことを示した。
(2) 励起状態の分子構造・断熱励起エネルギー

SAC‑CI 法は様々な性質の励起状態を精密に記述し、解析的エネルギー微分法が利用できること
から、励起状態の平衡構造や断熱励起エネルギーを詳細に研究することができる。HAB 型の分子を系統的に研究した例を示す。表１に HCO 分子、表２に HSiX(X=F, Cl)分子の結果を示す。HCO は開殻系分子であり、励起状態は２電子過程で記述されるが、general‑R 法は励起に伴う構造変化や断熱エネルギーをよく再現している。HSiF と HSiCl の１電子過程の励起状態は、SD‑R 法で良好な結果が得られる。励起に伴う分子構造の変化は静電力理論で定性的に予測することができるが、
SAC‑CI 法によるこれらの HAB 型分子の結果は全て静電力理論と合致するものである。
表１．HCO
表１．HCO の基底状態・励起状態の構造と断熱励起
エネルギー
状態

方法

励起
レベル

X2A' general‑R

1

Exptl.
A2Π general‑R

2

Exptl.
2

B A' general‑R
Exptl.

2

RHC

RCO

(A)

(A)

∠HCO

表２． HSiX (X=F, Cl)のCl)の基底状態・励起状態の構造と断熱励起エネルギー
分子

状態

方法

励起
レベル

RSiX
(A)

RHSi
(A)

X1A'

SD‑R
Exptl.

0

1.615
1.605

1.522
1.530

96.9
97.0

‑
‑

A1A"

SD‑R
Exptl.

1

1.613
1.609

1.541
1.484

110.5
111.0

3.113
2.884

X1A'

SD‑R
Exptl.

0

2.080
2.064

1.514
1.561

96.2
102.8

‑
‑

A1A"

SD‑R
Exptl.

1

2.058
2.047

1.516
1.499

113.6
116.1

2.707
2.569

Te

(degree) (eV)1.116 1.174

126.3

1.110 1.171

127.4

1.061 1.183

180.0

0.983

1.064 1.186

180.0

1.153

HSiF

0.000

1.112 1.368

105.9

4.741

1.160 1.362

111.0

4.797

∠HSiX
Te
(degree) (eV)HSiCl(3) 内殻イオン化スペクトル
内殻イオン化状態は多くの研究が報告されているが、軌道緩和が大きく、理論的に様々な性質の状態を一様に精度良く記述することは必ずしも容易ではない。とくに、内殻イオン化で観測されるサテライトピークは、軌道緩和と電子相関の効果が共に大きい。C1s の内殻イオン化エネル
ギーを系統的に研究した結果を図２に示す。SAC‑CI 法は実験で観測された化学シフトを精密に再現していることがわかる。また、図３には CH4 分子の C1s イオン化状態のサテライトスペクトルを比較する。理論は実験をピーク位置および相対強度ともによく再現していることがわかる。理論計算から、ピーク 1,2 は 1t2‑3p,ピーク 3 は 1t2‑3d, 4d, ピーク 4,5 は 2a1‑3s, 4s のリドベルグ励起が関与する shake‑up 状態であることが示された。
298

CO

2

2

(a) SAC-CI general-R

297

CH4, C1s

CO
296

CH F

295
CH Cl
2

2

CH F
3
HCN

294
293

H CO
2

HCONH

4

Intensity(a.u.)

SAC ‑CI (eV)

2 2

2

CH Cl
3

1

1
CH

5
(b) ESCA

2

292
291

3

3

4

5

C H

2 2

4

C H

2 4

290
290

291

292

293

294

295

296

297

290

298

300

Exptl. (eV)

図 2. C1s の内殻イオン化エネルギー：
SAC‑
SAC ‑ CI と実験の比較(4) 内殻イオン化状態における構造緩和
内殻イオン化状態における構造緩和
内殻イオン化状態において分子は構造緩和する。SAC‑CI general‑R 法の構造最適化によ340

表 3. 内殻イオン化状態の構造と断熱イオン化ポテ
ンシャル
分子
状態
CH4
基底状態
C1s イオン化
状態

計算し、実験結果と比較した。CH4 および NH3
NH3

基底状態

化により、C‑H 結合距離が 0.05A 短くなること
が実験の振動スペクトルから示唆されている

330

図 3. CH4 分子の C1s サテライトスペクトル(a) SAC‑
SAC ‑ CI と(b) ESCA の比較

り、平衡構造や断熱イオン化ポテンシャルをの結果を表３に示す。CH4 分子では内殻イオン310
320
Energy(eV)N1s イオン化
状態

方法
SAC‑CI
Exptl.
SAC‑CI
Exptl.

Re(A)
1.08441.08701.03401.0390

∠HXH
Td
Td
Td
Td

SAC‑CI
Exptl.
SAC‑CI
Exptl.

1.0082
1.01380.9820107.27107.23118.75IP (eV)

290.50290.86405.15405.52

が、本結果はこれを良く再現した。NH3 分子では、ほぼ平面構造に緩和し、N‑H 距離は CH4 分子と同様に収縮する。これらの分子の構造緩和は、
静電力理論により説明することができる。また、H2O 分子では O1s イオン化状態のポテンシャル曲面を計算することにより、振動スペクトルを精密に再現することに成功した[4]。
謝辞 本研究は、中辻 博 教授（京都大学）、倉本 圭 君(京都大学) 、石田真弓 博士(大正製薬)、
豊田和男 博士(大阪市立大学)との共同研究により行なわれました。また、内殻イオン化の研究においては、上田 潔 教授（東北大学）との共同研究により理論と実験の直接比較の研究を実現できました。深く感謝致します。
[1] H. Nakatsuji, Chem. Phys. Lett. 177, 331 (1991). [2] M. Ehara, M. Ishida, K. Toyota,
H. Nakatsuji, Reviews in Modern Chemistry, p. 293‑319 edited by K. D. Sen (World Scientific,
Singapore, 2002). [3] M. Ehara, S. Yasuda, and H. Nakatsuji, Z. Phys. Chem. 217, 161‑176(2003). [4] R. Sankari, M. Ehara, H. Nakatsuji, Y. Senba, K. Hosokawa, H. Yoshida, A. De
Fanis, Y. Tamenori, S. Aksela, K. Ueda, submitted.

2Tp05
Fragment MO 法によるタンパク質の丸ごと計算（産総研・計算科学）北浦和夫【序】コンピュータの高性能化に伴って、数千原子からなる分子・分子系の ab initio MO 計算が可能になってきた。ab initio MO 計算は、系のサイズの 3 から 4 乗に比例して計算資源（計算時間やメモリ）が必要になることから、系が大きくなると急激に計算が困難になるが、最近、密度汎関数法（DFT）
による 2000 原子系の計算が報告された。１） 一方、より高速に、より大規模な系を計算するために、
いわゆる linear scaling または order N 法と呼ばれる方法の開発が進められている。２，３）これらは計算資源が系のサイズに比例する程度でおさまる方法であり、これにより巨大系の電子状態計算が容易なることが期待されている。一方、巨大分子・分子系を計算するための近似法の開発も進められており
４）、我々も、フラグメント分子軌道法（FMO）を提案してきた。５）本発表では、FMO 法の現況と今後の開発計算について述べるとともに、FMO 法によるタンパク質の計算例をいくつか紹介する。
【FMO 法】FMO 法は、分子や分子集合体をいくつかの小さなフラグメントに分割し（図１）、フラグ
メント（モノマーと呼ぶ）とフラグメントペア（ダイマーと呼ぶ）について ab initio MO 計算を行うだけで、全系のエネルギーとプロパティが計算できる方法である。モノマーとダイマーの計算は、通常の ab initio MO 計算とほぼ同様である。Hartree-Fock（HF）レベルでは、次式を解く。

~~
F x C x = S x C x εx ,

(1)

~
~
Fx = Hx +Gx ,~xx= H µν
+ Vµνx + ∑ Bi µ hi hi ν ,
H µν

(2)

図１ 分子の分割

ix はモノマー(x=I)またはダイマー(x=IJ)を示す。射影演算子の項は、フラグメントの境界にある原子の基底関数を、各フラグメントに割り振る働きをする。また、V は次式で示す静電ポテンシャルである。

K(µν λσ )
Vµνx = ∑  ∑ µ (− Z A r − rA ) ν + ∑ Dλσ
K ( ≠ x )  A∈K
λσ ∈K


(3)

これらの計算で、通常の HF 法と異なるのは、例えば、モノマー I の計算には、I 以外のすべてのモノ
マーからの静電ポテンシャルが含まれることである。したがって、モノマーの計算は、全てのモノマ
ーが self-consistent になるまで繰り返し計算を行う。ダイマー IJ の計算には、IJ 以外の全てのモノマ
ーからの静電ポテンシャルを含めて解く。分子の全エネルギーは、モノマーとダイマーの全エネルギ
ー、Ex、

{ (

1
~
~
E x = Tr D x H x + F x
2)} + E ,NR
x(4)を用いて、次式により計算する。

E = ∑ E IJ − ( N − 2)∑ E I .
I >J

(5)

I

ただし、式(4)の ENR は核間反発エネルギーである。双極子モーメントなどプロパティも、モノマーとダイマーのプロパティを用いて式(5)と同様に計算できる。
実際の計算では、高速化のためにいくつかの近似を用いる。一つは、遠く離れたモノマーからの静電ポテンシャルを Mulliken の AO 電荷や原子の点電荷で近似する。もう一つは、十分に離れたダイマ
ーは HF 計算せずに、単に、モノマー間の静電相互エネルギーで近似する。これらにより、ほとんど精度を落とすことなく、高速計算が可能になる。
FMO 法は、一般的にフラグメントのサイズを小表１ FMO 法と ab initio MO 法の全エネルギーの比較（HF/STO-3G レベル）5)

さくすると誤差（通常の abmolecule

PDB

initio MO 計 算 結 果 と の

＃of
res

total energy / a.u.

error
kcal/molFMO

ab initio

差）が大きくなり、10~20

Alpha-1 (chain A)

3AL1

13

-4855.275923-4855.275472-0.3

原子程度（ポリペプチド

Alpha-1 (chain B)

3AL1

13

-4855.570969-4855.570053-0.6では、2 アミノ酸残基程

Alpha-1 (dimer)

3AL1

26

-9711.189715-9711.188371-0.8度）以上にとると、全エ[Pro , Leu ]crambin 1EJG

46-17779.503014-17779.502399-0.4

ネルギーのエラーは数[Ser22, Ile25]crambin

46-17777.148729-17777.148366

-0.2

22

25

1EJG

kcal/mol 程度となる。2 ア
ミノ酸残基を 1 フラグメントとして計算したポリペプチドの FMO-HF/STO-3G の全エネルギーは、
1kcal/mol 程度の誤差で ab initio MO 計算の結果を再現する。
（表 1）
FMO法は、高効率で並列計算

表２ができる。現時点で、ABINIT-MP

system（中野ら）およびGamess（Fedorov

Lysoら）に組み込まれており、両者ともMPIやDDI（GDDI）で並列化されている。前者はは既に公開されている6)。修正版Gamessによる、

Lyso
H2O
256
H2O
1024

FMO 法計算のタイミングデータ# of res# of atoms
128
2,036
same as
above
−
768
−
3,072

WF/Basis set# of basis
HF/STO-3G6,252
MP2/6-31G**
20,660
HF/6-31G*
3,328
HF/STO-3G
7,168

elapsed
time
36M
3D15H
15M
29M

CPU
# of CPUPIII (1GHz)128
PIV (2.2GHz)
24
PIII (1GHz)
128
same as above

タンパク質と水クラスターの計算
時間を表２に示す。これらはすべて1点計算の時間（経過時間）である。lyso（lysozyme-ligand複合体,2036原子）のHF/STO-3G計算は、128台のPCを用いると36分で計算できる。この計算の並列化効率は約80%である。
【FMO 法の現況と今後の計画】現在、ABINIT-MP では、HF 法のエネルギーとその勾配計算が可能である。また、ABINIT-MP を用いた MD 計算プログラム 7)も開発されている。修正版 Gamess では、
HF、DFT のエネルギーとその勾配にくわえて、MP2、MCSCF のエネルギー計算ができる。また、諸熊グループで開発された ONIOM 法様の、波動関数/基底関数の混合計算と IMOMM 計算も可能である。今後は、PCM などの溶媒モデルの組込みやいろんな波動関数での 1 次微分と２次微分の開発を進める予定である。
1) F.Sato,T.Yoshihiro,M.Era,H.Kashiwagi, Chem.Phys.Lett.,341 (2001) 645
2) G.E.Scuseria, J.Phys.Chem.A, 103 (1999) 4782
3). S.Goedecker, Rev. Mod. Phys.,78 (1999) 997
4) A.Imamura, Y.Aoki, K.Maekawa, J.Chem.Phys., 95 (1991) 5419.
5) T.Nakano et al.,Chem.Phys.Lett. 351 (2002) 475, K. Kitaura et al. , Chem. Phys. Lett. 312 (1999) 319
6) http://moldb.nihs.go.jp/abinitmp/
7) Y.Komeiji et al.,Chem.Phys.Lett., 372(2003)342

2Tp06

ProteinDF を基盤としたタンパク質量子化学計算システムの開発（東大生研計算科学技術連携研）○佐藤文俊、井原直樹、上野哲哉、大槻浩幸、
金城友之、田原才静、恒川直樹、西村康幸、牟田元、吉廣保、連宗旺（富士総研）西川宜孝、村松伸哉、岡田一俊（アドバンスソフト）稲葉亨、小池聡、長峰康雄、小池秀耀（東大生研計算科学技術連携研／アドバンスソフト）柏木浩【序】タンパク質は生物体の主要構成成分で、多数のアミノ酸が連結した巨大分子である。全てのタ
ンパク質の約 1/3 は何らかの形で金属イオンを持っており、構造の安定性へ寄与したり活性中心を形成したりする。また、タンパク質は生体内温度で反応が進行する精密分子である。このように、タンパク質の機能の研究には、巨大なペプチド鎖や活性中心の機構を精度よく調べることのできる解析手段が必須である。
実験的研究においては、種々の解析装置の進歩やその解像度の格段の向上、タンパク質精製技術の発展、遺伝子操作による機能解析法の登場などにより、大成功を収めている。一方、理論的研究では、これまで分子のサイズがネックとなり古典論やモデルによる解析方法が主であった。
しかし、タンパク質が関与する様々な反応は通常の分子同様、電子状態が重要な役割を担っている。そのため、今後の計算機の能力向上を見越して、タンパク質全体を考慮に入れた量子化学に基づく解析手段を整備する必要がある。
当グループでは、アミノ酸と金属イオンで同程度の定量性を持ち、タンパク質全体を分子として取り扱うことができる量子化学計算手法として密度汎関数法に注目し、プログラム ProteinDFを開発してきた。これを用いて 104 残基の金属タンパク質シトクロム c の全電子波動関数計算に初めて成功した実績がある[1]。現在、この ProteinDF をベースにして、タンパク質の精密なシミュレーションが実行できる量子化学計算システムを開発している[2]。ポストゲノム時代のタンパ
ク質の基礎ならびに応用研究に必要な機能を追加・拡張し、基盤を整備し、これを汎用レベルにまで完成させることが目標である。
本発表では、現在開発中のシステムについて紹介する。
【ProteinDF を基盤としたタンパク質量子化学計算システム[2]】当グループは ProteinDF をベースに、1,000 残基規模の超大規模密度汎関数計算や、数十残基のタンパク質の精密な分子動力学計算が実行できるシステムの開発を目指している。これらのシ
ミュレーション達成をサポートする機能の開発にも力を注ぎ、計算結果のデータベースも構築している。
本システムは 5 つのサブシステムから構成されている（図 1）。
(1) プロテイン・エディタ：本システムの実行・制御・解析に渡る全機能を統括する環境。以下の、2～5 の項目を統合する。グラフィカルなエディタ機能、タンパク質構造に有用な大規模分子グラフィックス、およびそのグラフィカル・ユーザ・インターフェースで構成した。様々なタイプのシミュレーションにおいて、ユーザにストレスを与えないように、共通の操作性、統一性のある表現をもって統括する。

統括環境
(1)プロ
テイン・
エディタ

ProteinDFシステム(2)自動計算法

QCLO法の高速化と拡張(3)構造最適化・
ab initio MD

量子分子動力学法

遺伝子データベース

PDBデータベース(4)大規模タンパク質計算

大型計算機環境に対応
ネットワーク(5)タンパク質波動関数
データベース

系統的計算とデータベース化

図 1：ProteinDF を基盤としたタンパク質量子化学計算システム概念図(2)自動計算法：タンパク質のシミュレーションをほぼ自動的に達成できる方法として、擬カノ
ニカル局在化軌道（QCLO）法を発案し[3]、本システムに組み込んだ。QCLO 計算の効率良い並列処理による高速化や、QCLO の性質を利用した溶媒分子の取り扱い、（複数の）活性中心を持つ複雑な電子状態を計算する機能も開発している。
(3)構造最適化・ab initio MD：X 線構造解析や中性子散乱、多次元 NMR などの実験で得られた構造、およびホモロジーモデリング、古典的な分子動力学法などで決定されたタンパク質の立体構造をもとに、タンパク質全体あるいは活性中心近傍における構造の最適化を ProteinDF で行う方法を開発している。同様に、時間発展のシミュレーションも構築する。
(4)大規模タンパク質計算：これまでの研究から、100,000 軌道のカノニカル分子軌道計算が倍精度演算限界であると推定した。大型計算機を使用して、このような計算が実行できるプログラムの開発を目指している。これは valence double タイプの基底関数で 1,000 残基規模のタンパク質に相当し、機能性に富んだ多くのタンパク質が計算対象に含まれることになる。
(5)タンパク質波動関数データベース：タンパク質の代表的な構造単位や、あるタンパク質に対して一連の遺伝子工学的な操作がなされた構造を選出し、これらの全電子計算や構造最適化計算を本システムによって系統的に実行している。得られたデータは、収集してデータベースとして公開する予定である。
それぞれのサブシステムについての詳細は本討論会にて別途報告するので、ご参考頂きたい。
本研究は文部科学省 IT プログラム「戦略的基盤ソフトウェアの開発」において実施された。
【参考文献】[1] Sato, F. et al,, Int. J. Quant. Chem. 63 (1997) 245-256; Sato, F. et al, Chem. Phys. Lett. 341 (2001) 645-651.
[2] 佐藤, 生産研究, 55 (2003) 245-253.
[3] Kashiwagi, H. et al, Mol. Phys. 101 (2003) 81-86.