1Pp025‫ݗ‬次元アルゴリズムに基づくオリゴペプチドのコンフォメ−ションӕ析（立教大理１ 産総研２ NTT 物性基礎研３）○家入 寛子１ 常盤 広明１ ସ嶋 х兵２ 寺前 裕之３【序論】

生体内で様々な機能を発現するタンパクࡐは、２０数種のアミノ酸から構成される‫ݗ‬分子であるにもかかわらず、ある一定の安定した三次元構造を保ちながら存在していることが知られている。一般的に生体内でのタンパクࡐの挙動は外界（熱浴）との環境の中での振る舞いであるため、常に揺らいでおり、ネイティブ構造、フォールディング過程の構造（モルテングロビュール状態）のӕ析は‫ؼ‬年のタンパクࡐ化学の重要Ӏ題である。すなわち、‫ݏ‬素や受容体などのタンパクࡐのフォールディング過程を知ることにより、副作用のリスクを負うことのない薬剤の開発や、ミスフォールディングによって引き֬こされるЋ伝病の治療が可能となることが期待される。
現在、X 線ӕ析やレーザーӕ析などの実験的な手法を用いたタンパクࡐの立体構造ӕ析が行われているが、それらの方法は、主に結晶系や固体に対するものであり、フォールディン
グの経時的変化をୈ৞することは困難である。一方で、最‫ؼ‬の理論化学の発展によりタンパ
クࡐの結合水や周囲の水をも同時に取り扱い、タンパクࡐ系全体の構造変化を経時的にୈ随できる分子動力学(Molecular Dynamics)法が広く普及してきた。
しかしながら、タンパクࡐ分子のように分子内水素結合や疎水結合をもつ分子に対する MD‫ڐ‬算では、それらの弱い結合を正確に取り扱うことが不可欠であるため、タンパクࡐに対する種々の力場が提案されてきたが、それらの適用範囲には制限があった。
そこで本研究では、対象化合物の஢子状態を ab initio 法でӕき、エネルギー勾配法により決定された force field を用いて、‫ݗ‬次元アルゴリズムによって、原子運動間のカップリングを考慮した MD ‫ڐ‬算（１）を行った。運動方程式のӕ法にはヴェルレ法を用いた。‫ڐ‬算対象はタンパクࡐのプロトタイプサンプルとして、タンパクࡐを構成するアミノ酸の中でもっとも簡単な Glycine および Alanine の５量体を取り上げた。
【方法】

Glycine および Alanine の５量体の初期構造は、それぞれ Glycine のα−helix お

よび Alanine のβ−sheet を基に、HF/3-21G レベル、‫ڐ‬算プログラム Gaussian98 を用いて決定した。実際の MD ‫ڐ‬算は、1step=1.0a.u.として実行した。‫ڐ‬算プログラムには、‫ݗ‬次元
アルゴリズムに基づく MD 法を組み込んだ並列化 HONDO5 を用いた。今回使用した‫ݗ‬次元
アルゴリズムでは、運動に混合性を持たせることで、Dynamics ‫ڐ‬算を効率よく実行することができる。

【結果】Fig.1 に 5000~8000step までの Glycine５量体のトラジェクトリーを示した。エ

ネルギーの変化は約 0.06a.u.の幅があり、分子内水素結合交替をくりඉしながら、より安定なβ-sheet form へと変化している。

energy(a.u.)-1103.74-1103.76-1103.78-1103.80-1103.82-1103.84
5000

5500

6000

6500

7000

7500

8000

steps
Fig.1 Energy alternation of Glycine pentamer

Fig.2 には、Glycine５量体の特徴的な３つの構造を示した。今回の‫ڐ‬算では、様々なエネ
ルギー変化に伴い、α-helix,β-sheet およびγ-turn に対応するいくつかの重要な構造を見い出すことができた。また Fig.3 は、Glycine５量体の中性分子に加えて、そのカチオンおよびアニオン状態に対するコンフォーマ−の 12000step(a.u)までのエネルギー分布をプロットしたものである。ほぼ均等な分布になっているが、分布がわずかにエネルギーの‫ݗ‬い右側にシフトしている。特にアニオン状態では、中性分子の場合とは異なり、エネルギー分布が左側にシフトしていることから安定したコンフォーマ−が多く存在することが示唆された。

一方、カチオン状態は、中性分子と同様な分布となった。エネルギー的には中性分子が最も安定で、ついで、アニオン、カチオンの順となった。このように Glycine５量体において、その஢子状態の違いによって、異なる挙動をとることが認められた。Alanine５量体の結果については当日報告する。

the number of conformers

7000

α-helix form

6000
5000

Neutral

4000

Cation
3000

Anion

2000
1000
0
1

2

3

4

5

6

7

relative energy

β-sheet form

γ-turnFig.2 Conformers of Glycine pentamer

Fig.3 Energy distribution of Glycine pentamer【参考文献】（１）大田原一成、下川信祐、寺前裕之

２００２分子構造総合討論会（神戸）1Pp026
DFT-MD によるメタンハイドレートのゲスト分子運動シミュレーション（金沢大理）○井田朋智、松本大輔、加藤信彦、水野元博、遠藤一央【序】メタンハイドレートはロシア、北米、また日本近海においても多量の存在が確認されている。その全量は熱量換算で石油に匹敵し、将来の代替エネルギー源として注目されている。しかし、常温・常圧では非常に不安定な物質なため、資源開発には分子レベルでの制御技術が不可欠である。現在、理論的研究において、古典分子動力学（MD）によるハイドレートの物性解析が主流となっているが、この方法では水やメタン分子の電子状態に起因する吸着・解離過程や分子内振動および回転運動を同時に計算することは難しい。また、二粒子間のポテンシャルは経験パ
ラメータに依存するため、混合ガスや新規物質に対するハイドレート研究も困難となる。
1985 年 Car-Parrinello[1]によって電子状態を考慮した分子動力学法（第一原理分子動力学）が提唱された。この方法は密度汎関数法（DFT）によって粒子間相互作用を考慮するため、分子間・分子内運動を同時に計算できる。しかし、多くの Car-Parrinello 法は平面波基底を用いるため、
メタンハイドレートの様な結晶単位構造が大きなものには適さない。本研究の目的はハイドレー
ト内のゲスト分子運動に注目し、DFT を用いた非 Car-Parrinello 型の MD 計算によりメタンハ
イドレート構造の安定性を解析する。
【計算方法】
・交換・相関ポテンシャルの検討

MP2
BLYP
PW91
PBE
LDA

分子間相互作用の評価に用いる DFT 計算

5検討した。孤立波基底による MD 計算では
基底関数重なり誤差（BSSE）を考慮しない
ために、大きな基底関数を用いることで過剰な安定化エネルギーを生じる。従って基底関数と交換・相関ポテンシャルの組み合わせにEnergy (kJ/mol)において、最適な交換・相関ポテンシャルを

0-5より、過剰安定化を打ち消す必要がある。
Fig.1 に様々な交換・相関ポテンシャルでのメタン－水分子間相互作用エネルギー曲線を示す。横軸の C-O 距離は、水分子の H-O-H
平面上に O…H-C を直線に配置したときの距離を現している。全ての DFT 計算は ADF
プ ロ グ ラ ム [2] で 行 い 、 基 底 関 数 と し て

-10

3

4

5

C-O distance (Å)
Fig.1 Intermolecular potential energy surfaces
between water and methane.

double-ζ型基底+分極関数（DZP）を用いた。
比較として MP2/6-311++G(2d,2p)に BSSE を考慮したエネルギー曲線も示す。最も MP2 のエネ
ルギー曲線と近い形を示した交換・相関ポテンシャルは BLYP であった。その他のポテンシャル
は全て過剰な安定化エネルギーを示した。他の水－メタン分子配置においても同様の結果を得た。
このことから DFT-MD 計算において、交換・相関ポテンシャルは BLYP を使用することとした。

・分子間・分子内運動のシミュレーション
DFT-MD 計算では分子運動のみ扱う古典 MDとは異なり、分子内振動も取り扱うために 1MD(a)

全体の時間発展を長くする場合、短い 1 ステッ
プ時間は非常に多くのステップを要請する。そこで、水・メタン 1 分子のみの運動シミュレー
ションを行い、最適な 1 ステップ時間を検討した。MD 計算には速度ベルレ法を用い、系全体

C-O distance (Å)

ステップの時間を短くとる必要がある。しかし、

43の初期運動エネルギーとして 5kcal/mol を与え、水とメタン分子の衝突過程をシミュレートした。

1.2

Fig.2 に水－メタン分子の CO 距離とメタンのの時間変化はどのΔt においても全く同じ軌跡となった（Fig.2(a)）。メタン分子の C-H 結合距
離変化は、どれも平衡結合距離（1.096Å）を中心に振動運動することが再現できた。しかし、
Δt ＝1.5fs において、水とメタンが衝突した後の振動運動が破綻した（Fig.2(b)）。さらにΔt が0.5fs 以下の計算も行ったところ、0.5fs と全く同

1
CH bond length (Å)

まで変化させ計算した。水－メタン分子間距離

1.2

1
1.2

晶構造（sI 型）中の水分子が作る籠型クラスタ
ーを用いた[3]。籠型クラスターは、メタン分子∆t=1.5fs

1.1
1

計算において、Δt は 1.0fs 以下とした。

モデル構造として、すでに報告されている結∆t=1.0fs

1.1

じ軌跡となった。これらのことから今回の MD

・ゲスト分子運動シミュレーション∆t=0.5fs

1.1

C-H 結合距離の時間変化を示す。シミュレーシ
ョンは 1MD ステップ時間（Δt）を 0.5～1.5fs

(b)

0

50
100
Time (fs)Fig.2 (a) Time dependence of the distance
between water and methane molecules,(b) Time dependence of the C-H bondlength in methane molecule.

1 つに対して水分子 20 個で形成されている。分
子動力学シミュレーションは速度ベルレ法を用い、粒子間相互作用はスレーター型孤立波基底（DZP）を使用し、交換・相関ポテンシャルは BLYP として DFT 計算から求めた。DFT 計算は
全て ADF プログラムを用いた。詳細な結果は当日報告する。
【参考文献】
1.R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
2.Amsterdam Density Functional (ADF); http://www.scm.com.
3.E.D. Sloan, Jr., “Clathrate Hydrates of Natural Gases, Second Edition, Revised andExpanded”, Marcel Dekker, Inc.: New York, 1998.

1Pp027
DFT-GB 法によるクライゼン転位反応系の溶媒効果（筑波大化）○髙田 真吾・守橋

健二【序論】 コリスミ酸(1)からプレフェン酸(2)へのクライゼン転位反応は、コリスミ酸ムタ－ゼによって反応が促進される。この反応は、生体内でフェニルアラニンやチロシンなどを合成するシ
COOH
O
HO

HOOC

CH2
C

H2
C

COOH

H

HO

キミ酸経路で重要な反応である。コリス
C
O

COOH

ミ酸からプレフェン酸への非酵素反応の活性化障壁を求める実験は、Andrews ら
によって行われており、水溶液中で 20.7

H2

1kcal/mol と報告されている[1]。また、密度 汎 関 数 法 (DFT) に よ る コ リ ス ミ 酸 の

COO
O
HO

OOC

CH2
C

H2
C

COO

H

HO

C
O

diaxial と diequatorial の構造やプレフェン

COO酸の構造、それぞれの 2 つのカルボキシル基からプロトンが解離した状態(1΄, 2΄)

Hの構造は求められており、1΄→2΄反応の

2'1'

遷移状態計算が Houk らによって行われている[2,3]。また、1΄→2΄反応の溶液中の活性化障壁を求める QM/MM 計算が Jorgensen らによって行われている[4]。しかし、連続誘電体モデルを用いた 1→2 および 1΄→2΄反応の水溶液中の遷移状態を求める DFT 計算は行われていない。本研究の目的は、溶媒効果を拡張 Born 式で取り入れた密度汎関数計算(DFT-GB)法を用いて、1→2 および 1΄→2΄反応反応の水溶液中の非酵素反応における活性化障壁を求めることである。
【方法】

全ての計算において基底関数は 6-31G*を使用した。DFT-GB 法では、交換相関汎関数に B3LYP を使用し（これ以後、B3LYP-GB と表す）、連続体の誘電率を気相中で ε = 1、水溶液中
で ε = 80 として計算を行った。
まず、1→2 反応のモデル反応系として、1、2 それぞれから、1 つの水酸基と 2 つのカルボキシル基をそれぞれ水素に置換した 3→

H
O
H

H

CH2
C
H
O

4 反応ついて、B3LYP-GB/6-31G* 計
H

HH

H

Ts.c

H

O

3

H2
C

H
H

H

4

算を用いて気相中と水溶液中で活性
C
O

H

H

HH
Ts.b

化障壁と最適化構造を求めた。この反応では、遷移状態として椅子型(Ts.c)とボート型(Ts.b)が考えられるので、その両方の計算を行った。
次に 1→2 反応の気相中と水溶液中の構造最適化を B3LYP-GB 計算で行った。この反応は水溶液中では、2 つのカルボキシル基からプロトンが解離したジアニオン状態で反応が進むと考えられる。従って、1΄→2΄反応についても同様に、B3LYP-GB/6-31G*による気相中、水溶液中の構造最適化を行い、活性化障壁を求めた。

【結果】 3→4 反応について、B3LYP-GB 計算によって得られた相対エネルギー∆E、気相中と水溶液中のエネルギー差 ∆∆E は表 1 のようになった。表 1 より、遷移状態では気相および水溶液中の両方で Ts.b よりも Ts.c の方がエネルギーが低く、椅子型の遷移状態が有利であることが確認できた。また、溶媒和による安定化は反応物より遷移状態で大きくなることが分かった。これは反応物より遷移状態で電子構造はより分極し、溶媒和による安定化を受けやすいためである。
3→4 反応の結果より、Ts.b よりも Ts.c の方がエネルギーが低かったので 1→2 および 1΄→2΄反応の遷移状態は Ts.c のみ計算した。1→2 および 1΄→2΄反応の相対エネルギー∆E、気相中と水溶液
中のエネルギー差 ∆∆E は表 2 のようになった。表 2 より、1→2 および 1΄→2΄反応でも 3→4 反応の結果と同様に、溶媒和による活性化障壁の低下が反応物よりも遷移状態で大きいことが分かる。
特に、1΄→2΄反応では溶媒和による活性化障壁の低下が 28kcal/mol と大きい。これは遷移状態では−負電荷の大きい酸素原子をもつ 2 つの COO 基が中央付近に集まるため、分子内の静電的な反発が大きくなり気相中で不安定となる。しかし、水溶液中ではそれが溶媒分子（誘電体）によって緩和されるためである。Ts.c の気相と水溶液中の B3LYP-GB 最適化構造を比較すると、結合が切れる C-O 結合は、気相中の 2.162Åから水溶液中では 2.228Åに結合が伸び、新しく生成する C-C結合は、気相中の 2.498Åから水溶液中では 2.470Åに結合長が短くなっている。計算された水溶液中の活性化障壁は 1΄→2΄反応で 19.2kcal/mol と実験値 20.7kcal/mol と良い一致を示した。
表13→4 反応における相対エネルギー

表2

1→2 および 1΄→2΄反応の相対エネルギー∆E(kcal/mol)と∆∆E(kcal/mol)
HF/6-31G*∆E(kcal/mol)と∆∆E(kcal/mol)

B3LYP-GB/6-31G*∆E∆E

HF/6-31G*

∆∆E

ε=1

ε=1

ε = 80

3

0

0

0

-

Ts.c

46.5

25.9

15.7

Ts.b

49.3

27.8

4

-20.3

-15.7

B3LYP-GB/6-31G*

∆E

∆E

∆∆E

ε=1

ε=1

ε = 80

1

0

0

0

-

10.2

Ts.c

52.9

29.6

18.6

11.0

17.0

10.8

1΄

0

0

0

-

-19.2

3.5

Ts.c΄

69.7

47.3

(19.2)(28.1)∆∆E = ∆E(ε=1)-∆E(ε=80) (kcal/mol)

実験値（水溶液中の活性化障壁）：20.7 kcal/mol

参考文献[1] P.R.Andrews, G.D.Smith, and I.G.Young, Biochemistry, 12, 3492(1973).
[2] O.Wiest and K.N.Houk, J. Am. Chem. Soc., 117, 11628(1995).
[3] P.Kast, Y.B.Tewari, O.Wiest, D.Hilvert, K.N.Houk, and R.N.Goldberg, J. Phys. Chem. B, 101,
10976(1997).
[4] M.P.Repasky, C.R.W.Guimarães, J.Chandrasekhar, J.Tidro-Rives, and W.L.Jorgensen,
J. Am. Chem. Soc., 125, 6663(2003)

1Pp028 シスプラチン錯体およびそのアコ置換錯体における溶媒効果(京大院工)○小野紘平・佐藤啓文・榊茂好【序】シスプラチン、cis-PtCl2(NH3)2 1、(cis-DDP=cis-diamminedichloroplatinum) は抗癌活性を持つが、trans-DDP は抗癌活性を持たないことが知られている。細胞内では cis-DDP は Pt-Cl 結合が切れた状態が安定となり、このとき溶媒である水分子が代わりに配位した加水分解型錯体(Scheme 1 の 2 および 3)に変化
していると考えられている

2+

+

NH3

NH3
H2O
Cl

Pt

Cl

NH3

-

NH3

Cl

H2O
NH3

Pt

H2 O

-

Pt

H2 O

NH3

Cl

Cl

OH2

1

2

3

Scheme 1(Scheme 1)。加水分解型錯体はがん細胞の DNA のグアニン塩基に結合してG

DNA 合成を阻害する(Scheme 2)。これらの白金

G

H3N
Pt

+

G

H3N

G

錯体では水中や細胞内では気相中と電子状態および分子構造が大きく異なるはずである。本研

Scheme 2

究では、DDP およびそのアコ置換錯体の構造、電子状態等に対する溶媒効果について PCM 法および RISM 法を用いて理論的検討を行った。
【方法】構造最適化計算は DFT 法を用いて行い、交換相関積分項は B3LYP 関数で近似した。Pt、Cl の内殻電子は有効内殻ポテンシャルで置き換え、それらの価電子には split valence 型基底関数を用いた。その他については 6-31G*を用いた。溶媒効果は PCM 法および RISM 法で評価した。
【結果】Figere 1 に gas 中及び水中におΔEkcal/mol14

ける cis-trans 間のエネルギー差(∆E)を示した。gas 中では trans-DDP の方が12

cis-DDP より 12.7kcal/mol 安定であった。

10水中でも trans-DDP の方が cis-DDP よtransが
安定

gas
12.7
PCM

8

RISM

6

り安定であるが、エネルギー差は

4

12.0kcal/mol とわずかに減少した。モノ

2

5.6

RISM
gas

trans 体で差がほとんど無いが、水中で

RISM

PCM

-2
cisが
安定

1.0
0.1

0

アコ錯体の場合は gas 中では cis 体と

12.0

gas

-4

は cis 体の方が 2.4kcal/mol 安定だった。

-6ジアコ錯体になると gas 中でも cis 体の

-8

-2.3

-2.4
-4.5
PCM
DDP

モノアコ錯体-7.9
ジアコ錯体

方が 4.5kcal/mol 安定となった。水中では cis 体がさらに安定であり、その差はFigure 1. gas中での構造を用いたシス体・トランス体の相対安定性(ΔE＜0のときcisが安定)

7.9kcal/mol であった。このことから cis体の方が水中での安定化が大きいこと、および水置換が進むにつれて cis 体の安定性が増大される
ことがわかった。
cis-DDP および trans-DDP の水和によるエネルギー変化について調べた(Table 1)。PCM 法を用いた場合、水中では cis-DDP の電子構造の Table 1. 水和によるcis-DDPおよびtrans-DDPのエネルギー変化(kcal/mol)
ひずみ変化が大きく、これによってエネ
Aqueous solution
PCM
RISM

ルギー不安定化が起こる(∆Edis)。一方、

Gas

溶媒和による安定化エネルギー(∆Esolv)

Egas

ΔEdis

ΔEsolvEaq

ΔE' solv

E'aq

は cis-DDP の 方 が trans-DDP よ り

cis

0

11.1

-36.0

-24.9

1.5

1.5

18.8kcal/mol 大きかった。また、RISM 法

trans

-12.7

2.4

-17.2

-27.5

8.6

-4.1

に よ る 水 和 自 由 エ ネ ル ギ ー (∆E'solv) はΔ(cis-trans)

12.7

8.7

-18.8

2.6

-7.1

5.6

cis-DDP では 1.5kcal/mol、trans-DDP では
8.6kcal/mol であった。次にアコ錯体の水和によるエネルギー変化について調べた(Table 2、3)。その結果、溶媒和エネルギーは水置換の進行に従って大きくなることがわかった。これは Cl-の脱離によって正電荷を帯びるためと考えられる。また、DDP では cis 体のほうが溶媒和による影響を強く受けるが、アコ錯体では溶媒和エネルギーは cis 体と trans 体いずれも同程度であった。

Table 2. 水和によるモノアコ錯体のエネルギー変化(kcal/mol)
Gas

Table 3. 水和によるジアコ錯体のエネルギー変化(kcal/mol)

Aqueous solution
PCM

Gas

Aqueous solution

RISM

PCM

Egas

ΔEdis

ΔEsolv

Eaq

ΔE'solv

E'aq

Egas

ΔEdis ΔEsolv

RISM
Eaq

E'aq

ΔE'solv

cis

0

5.7

-69.5

-63.8

-29.0

-29.0

cis

0

0.7

-202.6-201.9-148.8-148.8

trans

-0.1

7.0

-68.3

-61.4

-29.9

-30.0

trans

4.5

5.1

-203.6-194.0-151.0-146.5Δ(cis-trans)

0.1

-1.3

-1.2

-2.4

0.9

1.0

Δ(cis-trans)

-4.5

-4.4

1.0

-7.9

2.2

-2.3

実際の水中での分子レベルの情報を得るために RISM 法を用いて cis-DDP および trans-DDP の各原子に関する水の動径分布関数を調べた(Figure 2)。Cl 原子に対しては、2Å付近にピークが見えることから水の H 原子が近くに位置し、弱い水素結合を形成していると考えられる。ほかの原子に関しても検討し、Figure 3 のような溶媒和構造をとることがわかった。さらにアコ錯体についても同様に各原子に関する水の動径分布関数を調べた。Pt 原子に配位する水分子の数を計算したとこ
ろ、DDP、モノアコ錯体、ジアコ錯体でそれぞれ 18.9 個、19.9 個、20.8 個となり、加水分解の進行に従って、水分子の配位数が大きくなることがわかった。
H

H
O

4.62

3.97

g(r) 2.5

Pt

2
1.5

系列1
O
系列2
H
H 3
系列

1
0.5

H

4.94

Cl

N
H

H

3.00

2.60

0
2

4

6

8

O

10
r(Å)

Figure 2. Cl原子に関する水の動径分布関数

H
O

2.47

0

4.04

3.67

H

3.60

H

H

Figure 3. 各原子に対する水分子の配位の様子(単位はÅ)

1Pp029
分子動力学法を用いたアズリン蛋白の主鎖局所ダイナミクスの理論的研究（a 金城大社, b 金沢大理）○齋藤大明 a, 杉山歩 b, 吉本高志 b, 長尾秀実 b, 櫻井武 b, 西川清 b.
【序】タンパク質等に代表される生体分子の機能はマクロな分子の空間的な配置やそのダイナミク
スと密接な関係がある。例えば、プロテイン分子内の大きな変位をもつ原子の運動はマクロ分子の大きな部分の変位を促し、それらの振動特性には強い特定の 2 次構造依存があることが知られている。またタンパク質内の局所構造における結合やダイナミクスは 2 次、3 次構造を安定させる重要な役割を担っている[1,2]。このようにタンパク質のより詳細な構造の決定やそのダイナミクスの理解は、生体分子の機能解明への重要な因子であり、より精度の良い実験や解析法が望まれている。これらの機構解明へのアプローチとして分子動力学（ＭＤ）法がある。
分子動力学シミュレーションはタンパク質とその周囲の溶媒の詳細な運動の情報を時系列に得ることが可能であり、タンパク質の構造、ダイナミクス、そしてその機能性の関係を得るための重要なツールとなってきている。
本研究でとり扱うアズリンはブルー銅タンパク質の一種で植物やバクテリア等のさまざまな環境に存在する物質である。このクラスの比較的小さな金属タンパク質は、その銅原子を含んだアミノ酸の局所ドメイン（アクティブサイト）の原子配置や電荷などの環境に支配的である
ことが知られており、この銅原子やその周辺の残基の同位体効果によって分光学的特性も大きく変化することが知られている[3]。また銅原子の電子移動がたやすく起こる事から、自然において電子移動の代理を起こすプロテインであるとされ、その機構解明の研究は生物化学分野のみならず物理化学領域においても注目を集めている[4]。
本研究ではこのアズリン蛋白質活性中心であるアクティブサイトやその主鎖の局所ダイナミ
クスの理論的研究を分子動力学法を用いて行い、銅原子やタンパク質主鎖の振動特性や緩和パ
ラメータ等を導出し、実験等と比較可能な物理量の予測を行う。
【計算】
アズリンは 128 個のアミノ酸残基がペプチド鎖を形成することで構成されるタンパク質であり、それは８個の独立したβストランドと 1 個のαへリックス、そして 1 個の銅イオン原子(Cu2+)によって構成される[図１]。ここで銅原子を含んだアクティブサイトは図 1 において、アズリンタンパク質の一端に配置される。また、このアクティブサイト周辺の様子を図 2 に示す。同図において、Cu 原子は周辺のシステイン残基(Cys112)のＳ原子と２つのヒスチジン残基(His46, His117)のＮ原子と強く結合をつくり、メチオニン残基(Met121)のＳ原子やグリシン
残 基 (Gly46)の O 原 子 と 弱 い 相 互 作 用 を も つ と さ れ て い る [4]。 計 算 に 用 い た 初 期 座 標 はBrookhaven Protein Data Bank のデータベース“1AZU”を使用した。これは Pseudomonasaeruginosa (Cu2+) のＸ線構造解析から得られた座標であり、その分解能は 2.7Åである[5]。
MD 計算は Amber7[6]を用い、系の水和化、構造最適化、平衡化の操作を周期境界条件のもとで行った。アクティブサイトにおける原子の電荷分布や結合パラメータは量子化学計算により導出したものを用い、溶媒は TIP3 ポテンシャルの剛体水分子モデルを用いて系の温度を 300Kとして計算を行った。平衡化の操作は始めに NPT アンサンブルで 0K から徐々に温度を上昇させ、系の密度や体積、温度が一定の状態となるまで計算を進め、系が十分に平衡状態に達したのを確認したのち、NVE アンサンブルで計算を行いて解析に用いるデータのサンプリングを行った。解析した物理量は NMR 実験で得られる緩和パラメータと比較・検討可能な物理量である Spectral density とした。これを MD 計算によって得られるタンパク質のトラジェクトリー
から相関関数を計算することで導出し、NMR 等の実験の予測を行う。また各残基のオーダー
パラメータ S2 や平均自乗変位等の計算も行い、タンパク質分子の全体・局所的な回転や振動特性の検証を行う。詳細は当日報告する。

図 1：アズリン蛋白質の分子構造

図 2：アズリンのアクティブサイト構造【参考文献】
[１]
[２]
[３]

[４]
[５]
[６]

[７]

H. Frauenfelder, F. Parak, R.D. Young, Annu. Rev. Biophys. Chem. 17 451 (1988).
T. Alber, S. Daopin, K. Wilson, J.A. Wozniak, S.P. Cook, B.W. Matthews, Nature 330,
41 (1987).
R.S. Czernuszewicz, B.C. Dave, J.P. Germanas, “Spectroscopic Method in
Bioinograric Chemistry (Eds. E. I. Solomon and K.O. Hodgson) ACS Symposium
Series 692, 220 (1998).
H.B. Gray, E. I. Solomon, “Copper Proteins”(John Wiley, New York, 1981) 3, 1-39.
E.T. Adman, L.H. Jensen, Isr. J. Chem. 21 8 (1981).
D.A. Pearlman, D.A. Case, J.W. Caldwell, W.S. Ross, T.E. Cheatham, III, S. DeBolt, D.
Ferguson, G. Seibel & P. Kollman. AMBER, a package of computer programs for applying
molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations
to simulate the structural and energetic properties of molecules. Comp. Phys. Commun. 91, 1-41
(1995).
S. Pfeiffer, D.Fushman, D. Cowburn, J. Am. Chem. Soc. 123 3021 (2001).

1Pp030

電子動力学法による電子励起状態および電子移動反応に関する理論的研究(東大院理・東大院総合†)○天野倫子、重田育照†、高塚和夫†【序】分子の性質を正しく理解するためには、核のダイナミクスばかりでなく電子のダイナミク
スを考慮に入れる必要がある。現在のところ電子と核を個々に扱うダイナミクスについての研究例は多いが、2 つのダイナミクスを同等に取り扱う研究例は少ない。しかし２つのダイナミクス
が大きく相互作用している系(例えば非断熱遷移が起こる系、つまり複数のポテンシャル面が交わる領域)では、電子相関を正しく取り込み、かつ核のダイナミクスとカップルさせながら電子ダイ
ナミクスを解く必要がある。そこで本研究は、時間依存配置間相互作用(CI)理論に基づいた電子動力学法を開発し、電子ダイナミクスの観点から電子励起状態や非断熱過程を明らかにすることを目的とする。まず、固定した核配置のもとで計算した電子波束より定常状態の基底・特定の励起状態のエネルギー・電子波動関数を求める手法を確立する。続いて、この電子波束とカップルさせた核の“古典”トラジェクトリーを計算し、同じ時間スケールでの電子と核のダイナミクスを実現した。具体的には、本手法を、電荷移動と結合の組換えが同時に起こり非断熱過程の重要性が指摘されている化学反応、B+ + H2 → BH+ + Hに適用し、電子と核の運動の相互作用について詳細な解析を行った。
電子波動関数を、核座標の関数である
CSF(configuration state function)と時間の関数である配置係数 C I (t ) との積の線形結合で展開し、CI行列を用いてC I (t ) についての時間依存Schrödinger

0.8

CSFのpopulation |CI (t)|2【電子波束計算と定常状態のエネルギー･波動関数】

CSFのpopulationの時間発展として表した結果を図１に示す。図１で、例えば|221001|は軌道占有数に| 220002 |0.4| 201021 || 201201 |

0.2
0.1
0

0

5

で運動していることがわかる。
求めた電子波束の自己相関関数から、当研究室で開発されたエネルギースクリーニング法[1]を用いて| 221001 |0.5| 220002 |0.4| 220020 |0.3| 212001 |0.2| 221001 |

0.1
0

0

5

10

20

15

Time[a.u.]
0.35| 221001 |

(c)

0.3
0.25

| 220101 |0.2| 221010 |0.15| 220110 |0.1| 212001 |

0.05
0

0

5

10

15

20

Time[a.u.]

エネルギースペクトルを求め、基底状態、特定の励起状態におけるエネルギー(Table1)と電子波動関数

20

(b)

0.6

CSFのpopulation |CI (t)|2

く、速い周期や遅い周期が混ざり合った複雑な周期

15

0.7

よる電子配置を表す。各電子配置は他のいくつかの状態とエネルギーをやりとりして、単一周期ではな

10Time[a.u.]
CSFのpopulation |CI (t)|2Table 2 の脚注参照)で電子波束を計算し、それを

0.5

0.3

+ H2 → BH+ + Hを想定し、3 原子が直線構造をと合において、固定した核配置(分子構造においては| 221001 |0.6方程式を解くことによって電子波束を計算した。B+る、(a) B+ + H2、(b) BH2+、(c) BH+ + Hの３つの場(CSFの電子配置)

(a)

0.7

図1

CSFのpopulationの時間発展

(Table2)を計算した。この方法で求めたエネルギーと波動関数は、従来のCI行列対角化による方法で求めた結果とよく一致する。本手法を用いると、任意の高エネルギーの固有値と固有関数を、直交条件を付加することなく計算することが可能である。
【電子ダイナミクスと核の古典トラジェクトリー】
次に、核の運動とカップルさせて電子波束を時間発展させる。核のダイナミクスを考慮した場
合、核座標の関数である CSF は時間依存変化することになり、非断熱相互作用として CSF の核
座標微分を加える必要がある。その結果、電子の時間依存 Schrödinger 方程式は、核配置固定の場合とは多少異なり､

i=

⎡ e
⎛
∂
∂
∂
i=
ΦJ −
Φ Φ
C I (t ) = ∑ ⎢ H IJ
(R(t )) − ∑ R k (t ) ⎜⎜ ΦI
⎜
⎢
∂
∂
R
Rk I J
2
∂t
⎝
J ⎣
k
k

⎞⎤
⎟⎟ ⎥ C (t )⎠⎟ ⎥⎦ Jとなる。一方、核の運動方程式は、電子とのカップリングの効果として電子波動関数の核座標依存項を加え、⎡ i= ⎛
∂
∂
d
∂H
Rk =
ΦJ −
Φ Φ
= Pk + ∑ C I* ⎢ − ⎜⎜ ΦI⎢⎣ 2 ⎝⎜
dt∂Pk∂Rk∂Rk I J
I ,J

⎞⎟ ⎤
⎟⎥ C
⎠⎟ ⎥⎦ J

⎡ ∂
∂ ⎛⎜
∂
∂
d
i=
∂H
e
(R) − ∑ Pl
ΦI
ΦJ −
Φ Φ
Pk = −H IJ= −∑ C I* ⎢⎜⎢⎣ ∂Rk2 ldt∂Rk ⎝⎜∂Rl∂Rk∂Rl I JI ,J⎞⎟ ⎤⎟⎟ ⎥ C J⎠ ⎥⎦となる。式中の赤で示した項が核を動かすことによって取り入れられる項である。これらの式を同じ時間スケールで解くことにより、電子と核の同等に扱うダイナミクスを計算した。本手法を、
B+ + H2 → BH+ + Hに適用し、電荷移動の初期過程における運動の相関を解析した。結果の詳細は当日報告する。

参考文献

[1]

K. Takatsuka and N. Hashimoto, J. Chem. Phys. 103, 6057 (1995).

Table 1. Comparison of energy calculated by the conventional method and the present method
The ground state
The 1st excited state
The second excited stateCI
Present
CI
Present
CI
Present+(a) B + H2
-25.44356 -25.44356-25.13621 -25.13598-25.13323 -25.13322(b) BH2+(c) BH+ + H-25.39195-25.35520-25.39194-25.35520-25.15308-25.19875-25.15292-25.19837-25.10167-25.16838-25.10161-25.16772

Table 2. Comparison of configuration coefficient calculated by the conventional method and the present method
The ground state
The first excited state
The second excited state
CSF Coefficient
CSF Coefficient
CSF CoefficientCI
Present
CI
Present
CI
Present
2220000.954
0.954
2210010.925
0.928
2210100.930
0.930(a) 202020-0.123-0.1232200020.183
0.176211011-0.109-0.109202200-0.123-0.123
201021-0.116-0.1142200110.130
0.1282220000.944
0.944
2210100.945
0.942
221001
0.882
0.879(b) 220020-0.138-0.1392120100.117
0.125212001-0.110-0.117220200-0.138-0.1392110110.069
0.074
220002-0.334-0.3412220000.865
0.863
2210010.641
0.708
220002-0.135-0.130(c) 221100-0.174-0.1782120010.167
0.169
2211000.698
0.729
220200-0.384-0.385220101-0.562-0.5562202000.473
0.448
The calculation method and basis set : CISDT/6-31G++
Structure [Å] (a) r BH = 2.40, r HH = 0.75 (b) r BH = 1.20, r HH = 0.75 (c) r BH = 1.20, r HH = 1.50

1Pp031

時間分解光電子分光による
基底状態プロトン移動の実時間観測(東大院総合・東大院理) ○新崎康樹・山崎健太郎・高塚和夫【序】水素結合を持つ系におけるプロトン移動反応はもっとも基礎的な化学反応の一つであり、近年ポンプ・プローブ分光などによる励起状態ダイナミクスの研究がさかんに行われている。本研
究では、生体分子系などでより重要な意義を持つ電子基底状態でのプロトン移動ダイナミクスに着目する。電子基底状態と励起状態とでプロトン位置が異なる構造がより安定であることがしばしばある。この場合 2 つのポテンシャル面の差がプロトン位置により異なるために、電子励起状
態に用意した振動波束から、電子基底状態のプロトン移動障壁に対してより高い振動エネルギーまたはより低いエネルギーに選択的にダンブすることができると期待できる。そのようなダンプにより得られた振動波束を、遅延時間をフェムト秒のオーダーでコントロールしてイオン化する
ダンプ・プローブ光電子分光を用いることによって、プロトン移動ダイナミクスを実時間で観測できる可能性がある。前回までの発表で電子移動反応の実時間ダイナミクス [1] を探るのに用いた量子波束動力学計算 [2] と同様の手法で、クロロマロンアルデヒドを例にしたモデルポテンシャルを使ってこのような観測の可能性を調べた。
【モデル】ダンプ・プローブ光電子分光に関連する電子基底状態・励起状態・イオン基底状態のそれぞれについて、移動するプロトンの (x,y) 座標を座標系として選び、2 つの安定構造の近傍に合わせた調和振動子ポテンシャル V11 (x, y)、V22 (x, y) をおき、その間のエネルギー差に依存する
カップリング2
V12 (x, y) = f e−a(V11 (x,y)−V22 (x,y))(1)を設定し、これらを要素とする行列 V(x, y) を対角化することでモデルポテンシャルとした (図1(a)、(b))。安定構造のまわりの調和振動子のパラメタは ab initio 計算によりもとめ、カップリ
ングのパラメタ f と a は得られたモデルポテンシャルの遷移状態が ab initio 計算で得られる遷移状態構造にできるだけ合うように選んだ。このようにポテンシャルを取ることで、プロトン移動の前後におけるイオン化遷移確率の変化を安定構造のまわりでの値と関連付けた形でとりいれることができる。

(b)
1.0

O

y (Å)

y
H

(c)

xO

Left-Side Population Ratio(a)0.750 eV0.649 eV
0.546 eV

0.5

H

0.0

Cl

0.06

dump=5.26 eV
5.34 eV
5.42 eV

0.04

0.02

0.00
H

-0.5

0.0

x (Å)

0.5

1.0

0

100

200

300

Time After Dump (fs)

図 1: (a) 2 次元座標の選び方。プロトン移動の遷移状態を原点に取っていて、プロトン移動前では x > 0、移動後は x < 0 となる。(b) 基底状態モデルポテンシャル曲面の等高線図。(c) 3 つの初期条件 (ダンプ光エネルギー) での、x < 0 となる (つまり、プロトン移動した) ポピュレーションの割合の時間変化。
【量子波束計算】全系を関連する 3 つの電子状態で展開すると、ダンプ光 Vge・プローブ光 Vgk によってカップルされた、振動の波動関数 χ(x, y, t) についての連立運動方程式を得る (フェムト秒

ダンプ・プローブ分光の短かい期間内での回転は無視した)。
ih∂i χg (x, y, t) = T̂N + Vg (x, y) χg (x, y, t) + Vge (x, y, t)χe (x, y, t) +∂t

Z

dkk 2 Vgk (x, y, t)χk (x, y, t)
(2)

h
i
∂
∗
i χe (x, y, t) = T̂N + Ve (x, y) χe (x, y, t) + Vge(x, y, t)χg (x, y, t)
∂t
h
i
∂
∗
i χk (x, y, t) = T̂N + Vion (x, y) χk (x, y, t) + Vgk(x, y, t)χg (x, y, t)
∂t

(3)
(4)

ここで T̂N は運動エネルギー演算子、添え字の g、e、ion はそれぞれ電子基底状態、励起状態、イオン状態を、k は光電子の波数を表わす。イオンの連続状態を k での積分で表現した。これを離散化した kj (j = 1, 2, . . . , Nk ) での求積で置きかえ、kNk と Nk を積分が収束するのに十分な大きさにとると、上記の式中のハミルトニアンは (Nk + 2) × (Nk + 2) の行列 (Nk ≈ 30) となる。これを演算子分割法 (運動エネルギー演算子は FFT、ポテンシャル非対角項は解析的な対角化) により細かい時間ステップごとに対角化しながら時間発展をさせる。
【結果と考察】ダンプ光を 3 通りに選んだときの、x < 0 における振動波束のポピュレーションの時間変化を図 1(c) に示す。この 3 通りのダンプ光により、図 1(b) に等高線のある 3 つのエネル
ギーの波束を生じる。ダンプ光のエネルギーが 5.26 eV の場合はプロトン移動の障壁よりも高エネルギーの波束、5.42 eV では障壁よりも低い波束が生じるように設定してある。高エネルギー波
束ではダンプされた後に、プロトンが左右の安定状態間を振動している。この運動をイオン化によりプローブする。
例として、プロトンが左側に移動した場合のイオン化遷移確率が右側にいる場合よりも十分に大きい場合を示す。図 2(a) に、プロトン移動をおこす高エネルギーの振動波束の場合に得られた光
電子エネルギー分布を示す。ポテンシャル曲面の左側へと広い領域を通過する波束の運動に合せて光電子エネルギー分布に時間変化が見られる。図 2(b) はプロトン移動をあまりおこさない低エネルギーの振動波束の場合で、こちらでは光電子エネルギー分布には大きな時間変化は見られず、
ポテンシャル曲面の右側に留まっている波束から来るシグナルが見られる。

2.0
er

En

ct
le

0.5

ro

n

1.0
oe

250
200
Dela
y Tim
e (fs) 150

ot

300

gy

gy
er

1.5

Ph

oe

le

ct

ro

n

1.0
0.5
ot

250
Delay 200
150
Time
(fs)

Ph

300

En

1.5

(e

(e

2.0

V)

(b)

V)

(a)

図 2: 光電子エネルギー分布。(a) 高エネルギー波束、(b) 低エネルギー波束。
本発表ではプロトン移動反応の簡単なモデルを使って、ダンプ・プローブ光電子分光による基底状態プロトン移動の実時間観測の可能性を探った。プロトン移動の前後でイオン化遷移確率が大きく変化するなど良い条件が揃えば観測を十分期待できる。

[1] Y. Arasaki, K. Takatsuka, K. Wang, and V. McKoy, Phys. Rev. Lett. 90, 248303 (2003).
[2] Y. Arasaki, K. Takatsuka, K. Wang, and V. McKoy, J. Chem. Phys. 112, 8871 (2000).

1Pp032( B e + n H 2) @ C 6 0 の 動 的 電 荷 揺 ら ぎ に お け る
電子状態変化の解析（東大院総合）○重田育照・高塚和夫[ 序 ] 通常の金属内包フラーレンでは、内包されている金属がフラーレンを構成している炭素に比べて比較的重いため、運動が緩やかなのに対し、軽い原子の場合には骨格と同じくらいの時間スケールで運動すると予想される。それにより、構造の揺らぎの効果が大きく、内部の運動形態が複雑化する事が期待される。本研究では、Be 内包
フラーレンに水素分子を導入した系でのBeの位置が系にどのような影響をおよぼすかを明らかにするため、C60 内での Be の電子状態の変化を調べる。さらに Be 原子の動的な振る舞いを調べるため、有効電荷の時間変化を半経験的QM/MM分子動力学法 により解析し、内包される水素分子数や温度と有効電荷や内部運動の関係を明らかにする。
[ 計 算 結 果 ] 図１に、Be@C60 の Be を中心から５員環、６員環それぞれの中心方向（z軸とする）に動かした場合のポテンシャルエネルギー曲線を示した（B3LYP/6-31G(2d,2p)基底）。ただし、フラーレンの骨格の対称性は Ih に定め、結合距離を RC=C=1.367Å、RC-C=1.453 Åに固定した。Be の位置は中心付近に極小を持ち、ポテンシャルは単調に増加するのに対し、Beの有効電荷は最初減少し、1.2Å付近で増加に転じる事がわかる。つまり、Be が中心より大きく外れて運動することで有効電荷が大きく変化することが示唆される。R=0.0 Åの時、Be は C60 内でほぼ孤立系での電子状態を保っている。実際、1s 軌道は 61 番目の分子軌道、2s 軌道は最高占有軌道（HOMO）に大きなポピュレーションを有している。図２には HOMO と孤立系で計算した Be の原子軌道との重なりを示した。電荷が増加しはじめた 1.2 Å以降、2s 軌道の寄与は大きく減少し、その代わりに最初は2pz 軌道が、ついで2px・2py 軌道の寄与が増加してくることがわかる。つまり、Be は fullerene から受ける有効電場によって分極していることが分かる。

0.1

Z (6,B3LYP)

2s
2pz
2px
2py

1.4
1.2

Be

(/e)

0.05

Be

0
-0.05

Population (arb.)

0.15

E(5,B3LYP)
E(6,B3LYP)
Z Be(5,B3LYP)Z∆E (/ Hartree)

0.2

1
0.8
0.6
0.4

-0.1
0.2
-0.15
0
-0.2

0

0.2 0.4

0.6

0.8

1

1.2

1.4

1.6

RBe-COM (/Å)

図 １ Be@C60 分子のポテンシャルエネルギー曲線（E）と Be の有効電荷（ZBe）の変化(B3LYP/6-31G(2d,2p)) ５および６は５
員環、６員環の中心方向を示す

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

Distance from center of C60 (Å)

図 2 Be@C60 内での Be の原子軌道の変化
1s、2s、2pはそれぞれ孤立Be 原子の原子
軌道と、C60 内での HOMO との重なりを
表している

0.4

0.4
1,2,3,4,5
6,11,16,21,26
7,9,12,14,17,19,22,24,27,29
8,10,13,15,18,20,23,25,28,30
31,32,33,34,35
36,41,46,51,56
37,39,42,44,47,49,52,54,57,59
38,40,43,45,48,50,53,55,58,60
Be

0.3
0.2
0.1
0
-0.1

0.3
0.25
0.2
0.15
0.1
0.05

-0.2
-0.3

n=1, T=100(K)n=1, T=200(K)n=1, T=300(K)n=2, T=100(K)n=2, T=200(K)n=2, T=300(K)0.35Distribution (arb.)

Effective charge on equivalent C atoms (/e)

一方、Be の位置を変えても 1s 軌道の寄与が殆ど変化しないことから、Be の電荷の揺らぎは 2s、2p 原子軌道の寄与の大きな HOMO の変化と、Be の軌道の他の軌道へのし
み出しが原因である事が判明した。図３には、Be を５員環方向に動かした場合の炭素原子に誘起される電荷の和（対称性により同じ位置関係の電荷を足しあわせたもの）を示した。Be が動くことによって、ある特定の５員環を構成する全ての C 原子（１〜５）が負に帯電するのに対し、他の C 原子はそれを打ち消すように正に帯電している様子が伺える。Be の電荷が増加し始める 1.2 Å前後では他の C 原子も減少し始めることから、Be の電荷揺らぎは系全体の変化に支配されていることが結論付けられる。
一方図４には、AM1 QM/MM MD 法により得られた電荷分布を示した。水素分子の数が１つの場合には、有効電荷は温度を変えてもほとんど変化せず、一定の値（約ZBe=0.22）を示している。一方、水素分子の数をもう１つ増やすと、系の内部は実質３体の運動となり、非常に複雑な運動となる。図から分かるように、n=2の場合には全ての温度領域において電荷揺らぎが非常に大きく、さらに温度を変化させることによって、その分布が著しく変化することが分かる。特に n=2 の場合、低温 100K に比べて高温
200 〜 300K では、電荷分布のピーク位置が負の方へ大きくシフトしている。この原因を調べるために、C60 内での Be の軌跡を解析した。n=1 の場合には、Be は 5 員環と６員環の区別なくC60 の内部を運動しているのに対し、n=2の場合ではBeは６員環に沿って運動している、つまり５員環の中心付近には近寄らないことが明かとなった。さらに、温度が低い場合には、ある６員環の近傍でのみ運動している、つまりローカルミ
ニマムにトラップされているが、温度を高くするとバリアーをこえて隣接する６員環へと遷移することが判明した。その際、中心付近を長時間滞在するため電荷は負に帯電しやすく、電荷分布は負に片寄ることが予想される。つまり、大きな電荷揺らぎの原因は Be と C60 の中心方向を結ぶ運動に起因し、温度変化によるピーク位置の負の方向へのシフトは隣の６員環への遷移によるものであることが判明した。

0

0.2 0.4

0.6

R

0.8

Be-COM1

1.2 1.4

1.6(/Å)

図 3 炭素原子上の有効電荷の和（5 員環方向）
同じ関係にある炭素原子の電荷を足してプロットしてある

0-0.3 -0.2 -0.1

0

0.1

0.2

0.3

0.4

0.5

Z
(/e)
Be
図 4 AM1 QM/MM MD 法で得られた有効電荷分布（∆t=0.2fs、2ps の軌道を捨て、8ps の分布をしめす。温度制御には Nosé-Hoover 法を用いた）

1Pp033

分子動力学法による温度とストレスの制御とその応用(原研・産総研) ○志賀基之・篠田渉・三上益弘・蕪木英雄
一般に、金属や高分子などの固体材料の強度は、外部から系にある大きさの応力（ストレス）が加わった場合に固体内の原子間の化学結合を保持することができなくなるような限界値（降伏応力という）で特徴づけられる。定常な外部応力と歪みの大きさを関係づけ
る応力－歪み曲線には二つの領域がある。一つは応力の小さい領域で、応力と歪みが比例し（フックの法則）可逆的な弾性変形を表している。もう一つは応力の大きい領域で、フックの法則を逸脱した非線型な変形を示す。この領域では可逆的な変形に加えて不可逆的な塑性変形が含まれる。たとえば金属の塑性変形を原子レベルの立場から見れば、結晶中の欠陥や転位の運動をともなった結晶面のすべりによる原子の「ドミノ倒し」とそれによる結晶全体の構造変化が重要な要素である。
このように、材料強度の評価は外部応力下における化学結合のありかたを基礎に展開されるものであり、分子動力学・モンテカルロ計算や第一原理計算などの理論的アプローチ
から新たな知見が得られることが期待されている。このことを背景に、本研究では、一定の応力 t と温度 T と原子数 N を外部パラメータとした統計力学的集合（NtT アンサンブル）の分子動力学シミュレーション法について考える。その土台となるのは、Parrinello とRahman によって提案されたエンタルピーH を一定とする NtH アンサンブルの分子動力学法である[1]。本研究で用いる新しい手法では、実験と直接比較できる定温アンサンブルのものに変更し、さらに計算効率と精度を高めるために最近発展した方法（時間反転対称な RESPA 積 分 法 [2] に 基 づ く 四 つ の 定 温 サ ン プ リ ン グ 手 法 ： Nose-Hoover 法 、
Nose-Hoover chain 法、massive Nose-Hoover chain 法[3]、hybrid Monte Carlo 法[4]）を取り入れた。この方法で取り扱えるのは、応力を準静的にかけた可逆的変形である。すなわち、定常な応力下で熱力学的に安定な状態を保っている変形である。応力に対して線型だけでなく非線型の歪みを示す領域でも適用可能であるが、応力変化による不可逆過程である塑性変形の寄与は除外されている。
材料強度の基本的な評価には応力－歪み曲線のほかに、自由エネルギーの歪みに対する二階微分量である弾性定数行列 C も役に立つ[1,5]。特に、応力下での行列 C（または C に対して変形にともなう補正を施した弾性剛性係数行列 B ）は、その固有値がすべて正であることが系の安定性に必要な条件である（ボルンの安定性条件）ことから、降伏点において不安定になる原因（強く歪んで破壊がおこる方向）を解析する場合には有用である。このような事情から、Monte Carlo 法や分子動力学法で応力－歪み曲線と弾性定数を計算する方法が必要とされており、実際これまでの論文で幾つかのアプローチが紹介されている[6]。また、弾性定数については簡単な系で試行計算した結果から方法論の比較検討がなされている[7]。しかしながら、これらの量を精度よく求めるには多大な計算コストを必要とすることが知られており、効率的なシミュレーション手法の開発が望まれてきた。（とりわけ弾性定数は結果が収束するまでに通常およそ 1,000,000 ステップ以上かかってしまう。）
今回の 4 つの新しい方法を Ni 結晶 （ embedded atom potential, 500 原子, 300 K, 応力
t=0 の場合および一軸引張応力 t=5 GPa の場合）について弾性率計算に適用することで、それぞれの方法の計算効率（収束性）を調べた。その結果、 massive Nose-Hoover chain
法がすべての場合で最も優れていて、従来の方法と同程度の計算効率を与える
Nose-Hoover 法に比べて計算効率が 5-10 倍程度改善されることを確認した（図 1）。当日
は、この方法の詳しい定式化と応力歪み曲線の結果を報告する予定である。
t=0 GPa

t=5 GPa

1.5

C11 / eV Å-3

C11 / eV Å-3

Nose-Hoover
Nose-Hoover chain
massive Nose-Hoover chain
Hybrid Monte Carlo (10MD steps/MC cycle)

1.4

Nose-Hoover
Nose-Hoover chain
massive Nose-Hoover chain
Hybrid Monte Carlo

1

0.9

1.3
0

2

4
MD step

68[×106]

0

2

4
MD step

68[×106]

図 1：4 つの異なる NtT 分子動力学計算による Ni 結晶の弾性定数 C11 の収束性
[1] M. Parrinello and A. Rahman, J. Appl. Phys. 52 7182 (1981); ibid J. Chem. Phys. 76, 2662 (1982). [2] M. E.
Tuckerman et al., J. Chem. Phys. 97, 1990(1992). [3] G. J. Martyna et al., Mol. Phys. 87, 1117 (1992).
[4] B.
Mehlig, D. W. Heermann, and B. M. Forrest, Phys. Rev. B 45, 679 (1992). [5] M. Born and K. Huang, Dynamical
Theory of Crystal Lattices (Clarendon, Oxford, 1956); J. Wang et al, Phys. Rev. B 52, 12627 (1995); Z. Zhou and B.
Joos, Phys. Rev. B 54, 3841 (1996). [6] P. J. Fay and J. R. Ray, Phys. Rev. A 46, 4645(1992); J. V. Lill and J. Q.
Broughton, Phys. Rev. B 46, 12068(1992); I. Souza and J. L. Martins, Phys. Rev. B 55, 8733(1997); E. Hernandez, J.
Chem. Phys. 115, 10282(2001).
[7] M. Sprik, R. W. Impey and M. L. Klein, Phys. Rev B 29, 4368 (1983); J. R.
Ray, Comput. Phys. Rep. 8, 109 (1988); M. Karimi et al., Model. Simul. Mater. Sci. Eng. 5, 337 (1997); Z. Zhou and
B. Joos, Model. Simul. Mater. Sci. Eng. 7, 383 (1999).

1Pp034

結晶中における孤立ミュオニウム中心のダイナミクスに関する理論的研究（北大院理）○丸田悟朗・武田定【序】正ミュオン（ μ+ ）は、プロトンの 1/9 の質量、スピン 1/2，正電荷 e をもつ粒子であり、絶縁体や半導体などの物質中にミュオンを打ち込むと、その一部は電子と結合して孤立水素原子様の点欠陥とみなせる孤立ミュオニウム中心（Mu 中心）を形成する。物質中のミュオニ
ウムの振る舞いは、 μＳＲ法（ミュオンスピン回転／緩和／共鳴法）を用いて実験的に調べることができ、その結果から不純物中心としての孤立水素原子の構造や運動、反応について窺うことができる。一方、その質量が水素原子の 1/9 と軽く、トンネル運動などの量子効果が水素原子よりも際立って現れることから、孤立 Mu 中心の特異なダイナミクスに焦点を合わせた研究も、実験、理論の双方から精力的に進められている。
結晶中の孤立 Mu 中心の振舞いを研究するためには、結晶中での Mu の占める位置、すなわち Mu サイトを知ることがまず始めに必要であると考えられる。しかし、ミュオン実験だけから Mu サイトを決定することは極めて困難であり、相補的な理論的研究が Mu サイトを特定するために不可欠である。以前に我々はホウ酸および六方晶窒化ホウ素に対して μSR 実
験を行い、これらの結晶中で孤立 Mu 中心が生成して、何らかの運動をしている事を観測した。また予備的な理論計算により、これらの結晶中の Mu サイトを推測したが、その運動モードの詳細までは分からなかった。そこで今回、ab initio 分子軌道法/密度汎関数法によりこれ
らの結晶中における Mu 中心の断熱ポテンシャルを求めて、結晶内における Mu 中心の運動
モードを明らかにすることを試みた。
【計算】ホウ酸結晶を 12 分子の H3BO3 からなるクラスターで近似し、その中に Mu を置いて、基底関数は STO-3G から 6-31++G**レベルで、 UHF および UB3LYP 法を用いてGaussian98 で計算した。このときミュオンは静止しているものとして、水素原子核で置き換えて計算した。また格子間の Mu が引き起こす格子緩和は考慮していない。六方晶窒化ホウ素についても、水素で終端したクラスター(B12N12H12)2 を用いて同様の計算を行った。
【ホウ酸結晶】ホウ酸結晶中では、B を中心に３個の OH が平面三方型に結合した B(OH)3分子が二次元的に水素結合を形成して層をつくる。この層内には、水素結合環による空隙が存在する。二次元層の積層により、この空隙は一次元的なカラムを形成する。ab 面に垂直に z 軸をとり、z の関数として計算したポテンシャルエネルギーを図１に示した。ここで Muの x,y 座標は最適化されている。UHF/3-21G（□）では分子面の直上と直下に安定サイトがあり、面をすり抜ける時と上下の面に移る時にエネルギー障壁がある。UHF/6-31G*（●）では面をすり抜ける時の障壁がなくなる。電子相関を考慮した UB3LYP/6-31G(d')（■）では再び障壁が現れるが、今度は上下の面に移るときの障壁が低くなっている。用いる手法によりポテンシャルの形は大きく異なるが、いずれにせよ障壁の高さは非常に低く、零点エネルギ

ーを考えると、計算結果は Mu が非局在化したブロッホ状態にあることを示している。しかし、実験から 30K よりも高温側で Mu は熱活性型の運動をしていることが分かっているので、この温度領域では、ホウ酸結晶の格子が歪んで、Mu はスモールポーラロン状態にあると考えられる。
【六方晶窒化ホウ素】一方、六方晶窒化ホウ素（BN）は極端に異方的な結晶であり、強い共有結合で硬く結合した原子がつくる原子面が、弱いファンデルワールス力によって積層した構造を持つ。このため、打ち込まれたミュオニウムは原子面に挟まれた空間を２次元的に拡
散すると考えられる。z 軸を原子面に垂直にとり、z 座標を二つの原子面から等距離の高さに固定して計算したポテンシャルエネルギーを図２に示した。UB3LYP/6-31G(d')では ontop
サイトおよび hollow サイトが安定点と計算された。しかし実験によると BN 中の Mu の超微細結合定数は室温では真空中の値 A0 に等しいのに対し、ontop サイトでは 0.8A0 と計算された。
この系でも 70K より高温側では熱活性型の運動が観測されているので、ホウ酸結晶と同様に格子が歪んで Mu がスモールポーラロン状態にあると考えると、hollow サイトで計算された
超微細結合定数が 1.0A0 であることから、スモールポーラロンは hollow サイトに局在していると考えられる。しかし、実験結果から示唆されている低温でのコヒーレントな量子拡散には
ontop サイトが関与している可能性がある。
【結論】自己捕獲された孤立 Mu 中心の 1 次元（ホウ酸結晶）、２次元（BN 結晶）拡散が理論計算により支持された。より定量的な議論をするためには、格子緩和を考慮した計算が必要である。また極低温におけるコヒーレントな量子拡散の可能性が示唆された。

layersMu

図２ BN結晶中のMuのポテンシャル

図１

ホウ酸結晶中のMuのポテンシャル

1Pp035

SN2 反応における立体効果の本質(北大院工) 田地川浩人■序論■ 一般に SN2 反応は、反応の遷移状態付近での置換基の反転（Walden inversion）を伴うため、置換基の種類により、その後の反応性が大きく影響を受ける。たとえば、ハロゲン交換 SN2
反応
F- + R3CCl → Cl- + R3CF

（I）の場合、置換基（R）が水素原子（R=H）の時に比べ、メチル基(R=CH3)の反転を伴う SN2 反応の場
合、反応速度が数桁低下する。一方、反応の活性化エネルギーは、水素原子およびメチル基ともに負の活性化エネルギーを持つため、反応性の違いは動的立体効果が反応（I）へ直接的に影響していることになると考えられる[scheme 1]。
本研究では、SN2 反応における立体効果の本質を解明するため、置換基 R を水素原子（R=H）およびメチル基（R=CH3）について、ダイレクト・アブイニシオ・トラジェクトリー法により、理論的に研究した。

■計算方法■ まず、反応の活性化エネルギーを求めるため、Static な ab-initio および DFT 計算を行った。計算レベルは、HF/3-21G(d)から、B3LYP/6-311+G(d,p)および、QCISD/6-311+G(d,p)
まで、数種の計算レベルにて、R=H および CH3 について、それぞれの活性化エネルギーを計算した。反応座標に沿う、各点では、振動数計算により、平行点および遷移状態であることを確認した。
ダイレクト・アブイニシオ・ダイナミックス計算[1-3]は、まず反応種 CR3Cl(ここで R=H およ
および F-C-Cl
び CH3)の構造を HF/3-21+G(d)レベルで最適化し、C-F 間距離を 7.0-8.2Å の範囲で、角を 130-180 度の範囲でランダムに発生させ、初期配置とした。次に、CR3Cl の内部温度を 10Kになるように振動エネルギーを与え、２体間の衝突エネルギーを 0-30 kcal/mol とし、その後のダイナミックスを追尾した。メチル基の数の効果を明らかにするため、RnH3-nCl (n=0-3)について、
トラジェクトリー計算を行った。また、比較のため、仮想的に水素原子の質量をメチル基と同じにした系（m=15ｍH）とした場合のトラジェクトリーも計算した。

Scheme 1

■結果・考察■ Static ab-initio 計
算の結果、R=H および CH3 についての活性化エネルギーは、メチル基数（n）の増加とともに、増加したが、全ての計算レベ
ルで負の値を示した。
２体間の衝突エネルギーをゼロから
30kcal/mol までの範囲で 0.5kcal/mol づつ増加し、その後のダイナミックスを追
R=CH3（n=3）におけるトラジェクトリ
ー計算結果の１例を図に示す。時間ゼロ
ででスタートした F-は、0.1ps 後に錯合体を生成し、0.13ps でサドルポイントを

A-1

Energy / kcak mol

尾した。衝突エネルギー36.3kcal/mol、20

通過した。その後、Cl-イオンは脱離した。

10
0
-10
-20
-30

エネルギーの経時変化に見られるスパイ

-40
0.0

ク状のカーブは、メチル基の C-H 振動の

10がゼロの反応（F- + CH3Cl）は、衝突エネルギーがゼロの場合も反応が効率良く

0.3

0.4

0.5

0.6

8

Distance / A

クトリー計算の結果、メチル基の数(n)0.2B

複雑なカップリングによるものである。
R=CH3（n=0-3）についての、トラジェ

0.1

r1
r2

6
4
2

進行するが、メチル基数が増加するにつれて、反応の進行に必要な衝突エネルギ

0
0.0

0.1

0.2

0.3

0.4

0.5

0.6

ーが増加した。
質量 15 を持つ仮想的な水素原子によがゼロの場合においても進行した。
また、
アブイニシオ計算の結果、n=0-3 では、全て活性化エネルギーが負であることよAngle / degree

るトラジェクトリーは、衝突エネルギー

120

C

100

り、上記の特異的な反応性は、動的な立

θ1

80

θ2
φ

60

体効果であると結論できる。
0.0

0.1

0.2

0.3

0.4

0.5

0.6

Time / ps

Figure 1

ダイレクト・アブイニシオ・トラジェクトリー

法で計算した反応性トラジェクトリーの一例■References■(1) H.Tachikawa, M. Igarashi, and T. Ishibashi, J. Phys. Chem. A, 106, 10977 (2002)(2) H. Tachikawa, J. Phys. Chem.A, 104, 497 (2000).
(3) H. Tachikawa, J. Phys. Chem. A, 105, 1260 (2001)

1Pp036
変分崩壊のない 4 成分相対論的基底関数（ 九大院理・名市大） ○渡辺祥弘・舘脇洋【序】理論計算によって重い原子を含む分子を解析するには，相対論的量子化学が必要となる。
近年の理論や計算機の発展により 4 成分相対論的分子計算プログラムが整備されはじめてきた。
またこれらのプログラムを用いて高精度の分子計算を行うための基底関数の開発がなされている。
しかしながら，この基底関数の開発には，『 変分崩壊』が大きな障害となっている。この障害は重い原子になるにつれて顕著になる。Talman の，相対論的全エネルギー変分の最小点は，大成分にとっては最小点，小成分にとっては最大点であるという指摘からも，原子の全エネルギーを最小に最適化する従来の基底関数の生成方法では，変分崩壊が生じてしまう。
そこで，本研究では「変分崩壊のない」基底関数の生成方法の開発を目的とし，いくつかの重い原子について様々な even-tempered 基底関数を適用して変分崩壊発生との関連について調べた。これまでに得られた情報と基底関数について報告する。1【計算方法】今回の計算には，Dirac–Fock–Roothaan 法による原子計算プログラム 2 を用いた。
原子核モデルには一様荷電球模型を採用し，原子半径 R は Visscher らと同じく³´
R = 2.03952714 A1/3 + 1.39058668 × 10−5 bohr
（ A: 質量数）とした。また，計算結果の比較のため，Visscher らの数値 Dirac–Fock（ NDF ）解を厳密な解とした。基底関数には，even-tempered ガウス型基底関数を用いた。

ζi = αβ i−1(i = 1, · · · , m)

ζ1 = α = ζmax ,ζm = αβ m−1 = ζmin

— 基底関数の決め方 —1. 各対称性ごとに軌道指数の最大値 ζmax と最小値 ζmin を決め，これらを固定し，展開項数m を徐々に大きくしていく。
2. 1. で決めた ζmax と ζmin にそれぞれ別々の因子を掛け合わせて，軌道指数の範囲を変えて，同様に展開項数 m を徐々に大きくしていく。
【結果と考察】Hg 原子について計算を行った。軌道指数の最大値 ζmax は，各対称性において一番内殻の軌道のみを考慮した水素原子モデルの最適化計算より決定した。また，最小値 ζmin は，
既存の相対論的基底関数 3 と非相対論的基底関数から最も小さな軌道指数を選択し，安全を期して，0.7 の因子を掛け合わせた。以上の方法より決定した軌道指数の値は以下の通りである。
s+
p−
p+
d−
d+
f−
f+

max. exp. ζmax
4.94E+07
4.99E+07
4.02E+06
1.64E+05
3.30E+04
1.96E+03
1.55E+03

min. exp. ζmin
2.94E-02
2.45E-01
2.45E-01
1.14E-01
1.00E-01
6.52E-01
6.31E-01

各 s+ , p− , p+ , d− , d+ , f− , f+ 軌道について，それぞれ，7∼62, 6∼62, 6∼62, 4∼36, 4∼36, 3∼36,

3∼36 項の間で，展開項数を徐々に増やしていった。この時，得られた全エネルギーと NDF 解との差 ∆TE の変化を次の図に示す。
（ 以下，全ての数値には原子単位を用いている。）

0.00160.00140.00120.0010∆TE

0.00080.00060.00040.00020-0.0002-0.0004

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

set

項数を増やしていく途中に，NDF 解より下回っている部分を確認できる。次に，ζmax と ζmin にそれぞれ因子 A と B を掛け合わせて軌道指数の範囲を変更した。それぞれの因子から (A : B) と表現する。展開項数が徐々に増えていくに従って，∆TE は次の様に変化した。
set
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

∆TE(2:1/2)
1.2673076
0.2593685
0.0456886
0.0079311
0.0022886
0.0008131
0.0003413
0.0002611
0.0002121
0.0001723
0.0001444
0.0001207
0.0000996
0.0000812
0.0000652
0.0000514
0.0000396
0.0000296
0.0000213
0.0000145∆TE(5:1/5)
5.0249902
1.3039220
0.2722432
0.0543463
0.0106235
0.0026710
0.0008781
0.0001930
0.0000653
0.0000457
0.0000263
0.0000163
0.0000119
0.0000082
0.0000052
0.0000029
0.0000013
0.0000000-0.0000007-0.0000012∆TE(5:1)
1.4524751
0.3147354
0.0589784
0.0095323
0.0026789
0.0006042
0.0001522
0.0000916
0.0000435
0.0000314
0.0000233
0.0000169
0.0000126
0.0000095
0.0000071
0.0000053
0.0000045
0.0000034
0.0000030
0.0000029∆TE(10:1)
2.6111980
0.5946111
0.1025078
0.0225893
0.0054898
0.0008668
0.0003606
0.0000826
0.0000413
0.0000283
0.0000164
0.0000132
0.0000099
0.0000080
0.0000066
0.0000057
0.0000057
0.0000050
0.0000050
0.0000051

以上の結果より，変分崩壊は，大きな軌道指数を的確に与える事で回避される事が期待できる。
今回の計算で得られた，因子 (5:1/5) の 27 番目の組は，universal set となりえる。希ガス原子と

Es 原子について，Malli の universal set による計算との比較を以下に示す。
Ne
Ar
Kr
Xe
Hg
Rn
No
Es
1

(5:1/5)-128.691931-528.683762-2788.860563-7446.894038-19648.849251-23602.005522-36740.691809-33952.691620∆TE(27)-6.00E-11-1.24E-09-1.21E-08
2.66E-07-1.18E-06-2.73E-06-2.89E-07
1.71E-06

Malli-128.691931-528.683762-2788.860556-7446.893979-19648.848931-23602.005191-36740.691000-33952.690925∆TE(Malli)
3.28E-08
4.91E-07
6.55E-06
6.00E-05
3.19E-04
3.29E-04
8.09E-04
6.97E-04

H. Tatewaki and Y. Watanabe: J. Comput. Chem. in press
O. Matsuoka and Y. Watanabe: Comput. Phys. Commun. 139, 218 (2001)
3
T. Koga, H. Tatewaki, O. Matsuoka : J. Chem. Phys. 117, 7813 (2002)

2

NDF-128.691931-528.683762-2788.860563-7446.894039-19648.849250-23602.005520-36740.691809-33952.691622

1Pp037
静電場の可視化による新しい化学的描像（ＮＥＣソフト、ＮＥＣ情報システムズ１、ＮＥＣ基礎研２）○山口直人・半田享１・村瀬匡・高田俊和２【序】大型分子の計算においては、その膨大な数値情報の故に、非経験的分子軌道法の計算結果を、分子グラフィックスにより可視化することは、その物理的・化学的意
味を解釈する上で、極めて重要である。これまで、静電ポテンシャルの可視化が広く行われてきたが、今回、静電場の新しい表示法と必要な静電ポテンシャルの１次及び２次微分を計算するプログラムを開発したので、その基本的な考え方と、具体的な応用における有効性について、報告する。この結果から、静電場の形状から化学反応機構に関する理解が得られることが判明しており、タンパク質とリガンドの相互作用に関しても、新しい視点が開けるのではないかと期待している。
【静電場トポロジー描画法】静電場に限らず、ベクトル場には流れのトポロジーが存在しており、その様子はベクトル場の大きさがゼロになる特別の点から描くことができることが知られている[1]。しかしながら、これらの手法を分子に適応した場合、異符号を持つ荷電間に存在する静電的相互作用は、全く表現できないことが明らかになった。この問題を克服するため、今回、静電場のＣＧ表示を目的として、静電場が持っている基本的構造を幾何学的に解析して得られるトポロジー情報に基づいて描画する方法を開発した[2]。この方法の特徴は、１）ベクトル場の微分がゼロになる
特異点から静電場の特徴を抽出する、２）３次元空間に分布する適切な関数を用意し、その関数値を変更することにより、空間を選択的に表示する、の２点である。この方法により、同符号をもつ電荷はもちろんのこと、従来の方法では表示できなかった異符合をもつ電荷の作り出す静電場の骨格の表示も、可能となった。更に、ベクトル場
を表示するには、通常矢印や流線などを用いるが、それらが重なって表示されるため、化学的な特徴を見失ってしまう危険性があったが，それも回避できることになり、静電場の特徴を端的に示すことが可能となっている。
本描画法には、静電場（静電ポテンシャルの 1 次微分）と静電場の勾配（静電ポテ
ンシャルの 2 次微分）が必要であるが、それらの積分パッケージを独自に開発した。
これまで、マリケンチャージによる電荷を静電場の計算に用いていたが、これにより正確な静電場の表示が可能になっている。
【アミノ酸錯体への応用】
まず、本表示法の有効性を示すため、Mg-glycine-glycinemethylester の加水分解機能の解明に適応した。この分子は、マグネシウムにグリシンとグリシンメチルエステ
ルが配位する構造からなり、メチル基を除くとマグネシウムを中心にして対称な形をしている。しかしながら、実験からグリシンメチルエステルのみが OH-と反応し、加水分解反応を起こすことが知られている[3]。計算から見積もられた電荷からは、反応部位となるグリシンメチルエステルの炭素と、それと対称な位置にあるグリシンのカルボニル炭素との間で大きな違いはなく、どちらの炭素も同等に反応することが予測される。しかしながら、両者の周囲の静電的環境は大きく異なっていることが、本可視化手法によって、示されている。即ち、グリシンの炭素は電気力線で覆われ、OHが接近できないことが判明した。マリケン電荷による静電場と、分子軌道計算に基づ
く静電場とを比較して、図に示す。何れも、グリシンの炭素が電気力線で覆われていることが分かるが、分子軌道計算による結果では、グリシンメチルエステルの炭素近傍の静電場が広く開けており、OH-を誘導するような形状をしていることが示されている。
【終わりに】
タンパク質などの生体高分子の振る舞いに、その静電的な特徴が大きく反映されていることは明らかである。現在、In Silico スクリーニングと呼ばれているが，薬の候補分子の効率的なスクリーニングに、コンピュータの活用が望まれている。このような要請を受けて、色々な手法が提案されているが、信頼性の観点から量子化学計算への期待は大きい。本手法の活用分野として、In Silico スクリーニングが考えられるので、今後、タンパク質に応用して、その有効性を確認していく予定である。

図 1 点電荷による描像(左)、MO 計算による描像(右)参考文献[1] J.Helman, Lhesselink. Computer., 22, (1989) 8,27-36.
[2] S. Handa et al. J. Visual. Comput. Animat., 12, (2001) 167-180[3] S.Ina, Y.Hatano, H.Kashiwagi. Bioimages., 1, (1993) 1,37-40

1Pp038
Gauss-FE 混合基底法における基底関数の検討(豊田中研) ○山川俊輔、兵頭志明【緒言】 不均一場での電子状態計算を目的として、Gauss-FE(Finite Element) 混合基底法による電子状態計算を検討している 1)。本手法は Gauss 関数を用いた LCAO 法と比較して、原子数が増えても基底数の増大が少ない、基底の自由度が高いなどの特徴を有しており、不均一静電場が存在する状況下での分子やクラスターの計算に適していると考えている。Gauss 関数は原子核近傍での波動関数の表現能力を補うために組み込んでいるが、FE 基底との適切な組合せについては明らかではなく、基底の選択と計算精度との関連を把握しておく必要がある。そこで、本報告では CO
等の低分子を対象として、各基底の組合せと計算精度との関係について検討を行った。
【方法】 FE 基底としては Lagrange 補間多項式 2)を用いた。１次元の Lagrange 補間多項式は式(1)、(2)のように表現される。

ψˆ ( x) = N 1 ( x)ψ 1 + N 2 ( x)ψ 2 + L + N n ( x)ψ nn

N i ( x) = ∏j =1( j ≠i )(x − x j )( xi − x j )

-(1)

(a)

(b)

-(2)

1,･･･, n: 各要素内の節点番号

xi , ψ i , N i : それぞれ各節点の座標、関数値、基底関数
n=2 は線形補間で、図 1-a のような節点配置となる（以下、FE-L1 と称する）。n=3 では各軸 3 つの基底

図 1 直方体要素の節点配置(a) 8 節点モデル(Lagrange 1 次補間）(b) 20 節点モデル(Lagrange 2 次補間）関数が存在するため 1 要素内に 27 個の節点が必要であるが、ここでは高次項を除外した要素モデルを用い
a2て、1 要素内の節点数は 20 とした（以下、FE-L2 と称する）。節点配置は図 1-b のようになる。この FE 基底に Gauss 関数を組合せて混合基底とした。2 原子分子の計算に用いた要素分割の様子を図 2 に示す。a1 の値は原子核間位置に応じて変更した。a1 の 2 倍が核間距離となる。a2 の値は 1.0 bohr で固定した。x, y, z

a2

a11

a1

図 2 要素分割図(図中●印は原子核位置を示す）

各軸の要素分割数は 16 として、中心から最外縁の節点までの各軸上の距離は 20 bohr で固定した。
電子間クーロンポテンシャルは Gauss-FE 混合基底(FE 基底: Lagrange 2 次補間多項式)を用いてポアソン方程式を解くことによって評価した。ポテンシャルの境界値は各要素内の電荷分布を 4 重極子モーメントまでの多項式に展開して計算した。電子間の交換相関ポテンシャルには B-LYP を用いた。固有値、固有ベクトルは CG 法 3)を用いて計算した。
【結果】 まず核間距離を実験値(2.1316 bohr)で固定して FE 基底と Gauss 基底の組合せと計算精度との関係を検討した。比較項目を総電子エネルギーEe と双極子モーメントμとして表 1 に示す。図2 に示した要素分割は非常に粗く、原子核間に節点を一つしか含まないものであり、このような要素分割では、特に線形補間では要素内の波動関数の変化を記述できない。確かに、FE-L1 では、組合せる Gauss 基底の表現能力が低下すると近似

表 1 計算精度の比較の度合いは急激に低下する。他方 FE-L2 では、Ee / hartreeμ / Debye

6-31G (1)-132.1358

0.754

6-31G (2)-134.5675

0.593

6-311G (1)-132.5910

0.993

6-311G (2)-134.9489

0.514

6-31G (1)-135.2045

0.451

6-31G (2)-135.8467

0.150(2)+ FE-L2 基底を用いて、平衡核間距離 RCO、結6-311G (1)-135.5070

0.331

合エネルギーΔE、振動数νCO を評価した。表 26-311G (2)-135.8948

0.161

基底関数

FE-L1 に比べて近似レベルの低下は小さい。FE-L1
より表現能力の高い FE-L2 が Gauss 基底の不足部分を補っていることを示している。本検討の範囲

FE-L1

内においては、6-311G (2)+ FE-L2 が 6-311+G(d)と同等の計算精度を有していた。自由度の高い FE
基底によって全体の調整がなされていた。次に表
1 において、最も高い計算精度を示した 6-311G

FE-L2

より B-LYP/6-311+G(d)と同等の近似レベルを有し

6-31G のみ

-135.7638

0.104

ていた。

6-311G のみ-135.8055

0.105

6-311+G(d) のみ-135.85450.176

なお、表 1 の結果は C と O それぞれに属する
Gauss 基底間の重なりがないと計算精度は維持

Gauss 基底名の後の番号は以下の基底に対応されないことを示しているが、原子核近傍の要(1): 原子核から 1.0 bohr 以内で減衰する Gauss 関数を使用

素分割を検討することにより、実際には異核間(2): 原子核から 2.0 bohr 以内で減衰する Gauss 関数を使用の Gauss 基底が重なりを持たない条件でも波動関数を表現することができる。図 3 は、ビニ
表 2 CO の計算結果

6-311G(2)B-LYP/+ FE-L2

6-311+G(d)RCO / Å

1.127

1.139

1.128

基底の寄与部は各原子周囲に限られていΔE / kcal･mol-1

291.0

259.6

256.2

ることを確認できる。

νCO / cm-1

2193

2100

2143

実験値 4)

ルアルコールを計算した例であり、Gauss【まとめ】 本検討で用いたような粗い要素分割であっても要素内の補間関数に高次の関数を用いれば、核近傍のみに Gauss 関数を重ねることにより計算精度を維持できることを確認した。

H

H

C

C

O

【文献】

H

H

(1) S. Yamakawa and S. Hyodo, R&D Rev. Toyota CRDL,
38 (2003) 26.
(2) 鷲津久一郎

2ψ FE

他共編, 有限要素法ハンドブック I,(培風館, 1981) 138-140.
(3) C. K. Gan, P. D. Haynes, and M. C. Payne, Comput.
Phys. Commun., 134 (2001) 33.
(4) B. G. Johnson, P. M. W. Gill, and J. A. Pople, J. Chem.
Phys., 98 (1993) 5612.

2
ψ Gaussψ Gaussψ FE22
ψ total= (ψ Gauss + ψ FE )図 3 ビニルアルコールの電子密度分布(10-5 の等値面）

1Pp039

有効ポテンシャル解析接続法による実時間相関関数の計算（ 1 科技団・2 奈良女大理） ○堀越 篤史 1,2・衣川健一 2

0.4

0.4

0.2

0.2

potentials

potentials

有限温度量子力学系における実時間相関関数の計算法の一つに、セントロイド 分子動力学法 (CMD) がある [1]。これは有効ポテンシャル面上の古典運動についての分子動力学計
算から相関関数を計算する近似法であり、量子液体の輸送係数の計算等に用いられているのであるが、半古典的な手法ゆえ、量子コヒーレンスが重要な系においてはうまく機能しないことが知られている [2]。
CMD で用いられる有効ポテンシャルは effective classical potential と呼ばれるもので、
これは Legendre 変換で定義される通常の有効ポテンシャル (standard effective potential)とは異なる性質を持つが、ゼロ温度極限で両者は一致する。例えば二重井戸型ポテンシャ
1 4
ル系 V (q) = − 12 q 2 + 10
q におけるそれらの温度依存性は図 1 、図 2 のようになり、低温極限（ β → 大）で確かに両者は一致する。standard effective potential はどの温度においても常に convex であるが 、effective classical potential は必ずしもそうではない。そしてeffective classical potential が convex になっていない時、CMD は良い近似ではなくなる。
何故なら convex でないような、非常に非線形性が強いポテンシャルについての分子動力学計算（平均操作）は、系が有限温度下で保持している量子性を破壊してしまい、相関関数が過度の減衰を示してしまうからである。これは ensemble dephasing と呼ばれている[2]。
（高温極限の場合は系の量子性はもとから破壊されているので CMD は exact になる。）

0.0
−0.2

classical
β=0.1
β=1
β=10
β=100

−0.4
−0.6
−2

−1

0

qc

0.0
−0.2

classical
β=0.1
β=1
β=10
β=100

−0.4
−0.6
1

2

図 1: effective classical potential Vβc (qc )

−2

−1

0

Q

1

2

図 2: standard effective potential Vβ (Q)さて有効作用形式に従えば 、standard effective potential からダ イレクトに松原 Green関数（虚時間相関関数）が求められ、さらに解析接続を用いることで実時間相関関数を計算することができる。われわれはこれを有効ポテンシャル解析接続法 (EPAC) と呼び 、有限温度量子ダ イナミクスの計算手法として確立することを目指している [3]。EPAC の枠
組では、例えば２点相関関数は次のように計算される。
ÃCβAC (t)

=

h̄
βh̄ωβ
coth
2mωβ
2

Ã

!

cos ωβ t − i

h̄
2mωβ

!

sin ωβ t + Q2min

ここで ωβ は standard effective potential Vβ (Q) の最小点（ Q = Qmin ）における有効振動数である。EPAC は CMD の場合ような古典論的な平均操作を用いないため、CMD が破綻するような系においても有効に機能すると期待される。
今回は二重井戸型ポテンシャル系を例に
imaginary time
CMD との比較を通して EPAC の有効性をpath integral
議論した。経路積分分子動力学法を用いて Z = R Dqe−SE /h̄βeffective classical potential を計算し 、それから Legendre 変換を経てstandard effective
PIMC/PIMD∨
potential を求め、これら二種類の有効ポテ
effective
standard
ンシャルから CMD 相関関数、EPAC 相関
Legendreclassical> effective関数をそれぞれ計算した。計算の流れは右
potential transform potential
図のようにまとめられる。
Vβc (qc )Vβ (Q)
これら二種類の相関関数をexact な相関関classical
analytic数と比較したところ、両者とも短時間の振eq. of motion
continuation∨∨
舞いはexactなものに良く一致したが 、時間CMD
EPACが増大するにつれて両者ともにexactな相関
CANcACCβ (t) ' Cβ (t)Cβ (t) ' Cβ (t)関数からずれていった。ただしその振舞い
は、CMD 相関関数が ensemble dephasingCβ (ω) = E(ω) CβCAN (ω)
による過度な減衰であるのに対し 、EPAC
相関関数の方は振動を続けており、二重井
real time quantum><
戸型ポテンシャル系特有のコヒーレントな
correlation function
振動を近似的に再現することができていた。
Cβ (t) = hq̂(t)q̂(0)iβ
このような量子コヒーレンスが重要な系において、EPAC が有効に機能することが確
図 3: 計算の流れ認された。
本研究は科学技術振興事業団計算科学技術活用型特定研究開発事業「量子多分子系
ダ イナミクス・シミュレーションの確立と応用」に基づいて遂行されたものである。

[1] J. Cao and G. A. Voth, J. Chem. Phys. 99, 10070 (1993)[2] S. Jang and G. A. Voth, J. Chem. Phys. 111, 2371 (1999)[3] A. Horikoshi and K. Kinugawa, J. Chem. Phys. (2003) 印刷中

1Pp040
B-spline 関数の 2 原子分子への応用(中京大教養)

齊藤史郎

【序論】

B-spline 関数は，自由度の大きい区分多項式で数値計算上の線形従属性もかなり小さいという利点があり，原子波動関数の計算に応用されている．2002 年分子構造総合討論会において，我々は原子種，軌道角
運動量を問わず同一の B-spline 関数を用いて，核電荷が 54 までの原子とその ±1 価イオンについて高精度の Hartree-Fock-Roothaan 波動関数を得ることができることを示した．つまり，原子波動関数の計算においては，B-spline 関数は universal basis set となり得ることを示した．
我々は原子波動関数の計算に引続き，2 原子分子波動関数の計算に B-spline 関数を応用することを試みた．
ここでは，2 原子分子の波動関数は回転楕円体座標 (ξ, η, ϕ) で記述し，(η,ϕ) 成分の関数を球面調和関数で展開することによって Schrödinger 方程式の簡単化を行った．これにより，ξ 成分の関数は 1 次元 B-spline−1関数で展開できることになり，積分の計算量を減らすことができ，また，r12の Neumann 展開を有限項で2+のいくつかの対称性の異なる軌道に打ち切ることが可能になった．今回は，この方法を用いて H+
2 , HeH

B-spline 関数を用いて universal basis set としての有効性を考察した．
【理論と計算方法】
回転楕円体座標での運動エネルギー演算子 T と核−電子ポテンシャル Vnucl は次のようになる．
 2 
 


 
∂
∂
1
1
2
∂
∂
2
2 ∂
+
T =− 2 2
−
1)
)
,
(ξ
+
(1
−
η
+
R (ξ − η 2 ) ∂ξ∂ξ∂η∂η
ξ 2 − 1 1 − η 2 ∂ϕ2
2{(ZA + ZB )ξ − (ZA − ZB )η} .
Vnucl =2
R(ξ − η 2 )

(1)
(2)

ここで，R は核間距離，ZA と ZB はそれぞれ焦点 A と B 上の核電荷，1 ≤ ξ < ∞, −1 ≤ η ≤ 1, 0 ≤ ϕ ≤ 2πである．分子軌道 ψm (ξ, η, ϕ) は球面調和関数

Ylm (η, ϕ) = (−1)

m+|m|
2



を用いて

ψm (ξ, η, ϕ) =2l + 1 (l − |m|)!
4π (l + |m|)!



2

|m|

Pl

(η)eimϕCilm Xi (ξ)Ylm (η, ϕ)

(3)

(4)

il

と展開する．この関数を使用すると運動エネルギ−演算子は

 

∂
2
∂
m2
2
− l(l + 1)
T =− 2 2(ξ − 1)+ 2
R (ξ − η 2 ) ∂ξ∂ξ
ξ −1と簡単になる．式 (1) から式 (5) への変換では




∂
∂2
1
2 ∂
Y m (η, ϕ) = −l(l + 1)Ylm (η, ϕ)(1 − η )+∂η∂η
1 − η 2 ∂ϕ2 l

(5)

(6)

の関係を用いた．このように (η,ϕ) 成分の関数を球面調和関数で展開すると運動エネルギ−項が簡単になるだけではなく，(η,ϕ) での積分は 3j 記号で表すことができるようになる．

X(ξ) 関数は K 次 B-spline 関数 {Bi,K (ξ)}i=1,...,M を用いてXi (ξ) = (ξ 2 − 1)µ Bi,K (ξ)(7)

と展開する．使用した M 個の K 次 B-spline 関数 {Bi,K (ξ)}i=1,...,M を 1 ≤ r ≤ Ξ 内で構成するために，

M + K 個の節点 {ti }i=1,···,M+K を次のように生成した．

1, 
i = 1, · · · , K




i−K
cosh Ξ
ti =,i = K + 1, · · · , M
M −K +1 Ξ,
i = M + 1, · · · , M + K(8)

式 (7) の µ は通常原点での振舞の考察から m 対称性の軌道では µ = m/2 にとるが，今回はそれに加えて，universal basis set としての有効性を調べるために対称性に関わりなく µ = 0 とする計算も実行した．
【計算結果】

H+
2 の R = 2.0 au での各分子軌道のエネルギ−を表 1 に示した．B-spline 関数のパラメ−タは M = 31,
Ξ = 50, K = 9 にとった．µ = m/2 の計算結果は 2 次元 B-spline 関数の結果によく一致しているが，µ = 01の計算結果では π 軌道で大きなエネルギ−差を生じている．これは (ξ 2 − 1) 2 の B-spline 関数による記述が貧弱であることを示している．また，π 軌道以外では µ = m/2 と µ = 0 の計算はよく一致している．
当日は他の計算結果，多電子系への応用についても議論する予定である．
表 1. R = 2.0au における H+
2 の軌道エネルギ− (au)
今回の計算 aµ = m/2µ=0−0.602 634 214 495−0.167 534 392 202
0.071 228 180 108
0.071 282 298 247
0.273 300 373 359
0.273 300 842 084
0.287 267 318 209
0.287 267 318 191
0.375 037 457 080
0.375 037 457 065
0.376 874 499 299
0.376 874 499 205
a 1 次元 B-spline 関数
分子軌道
1sσg
2pσu
2pπu
3dπg
3dδg
4f δu
4f φu

Ref.[2]2 次元 B-spline 関数−0.602 634 214 494 95−0.167 534 392 202 38
0.071 228 180 104 14
0.273 300 373 356 34
0.287 267 318 189 24
0.375 037 457 063 22
0.376 874 499 203 25

Ref.[3]
有限差分法−0.602 634 214 497−0.167 534 392 205
0.071 228 180 106
0.273 300 373 37
0.287 267 318 24
0.376 874 494【参考文献】[1] 菅野, 吉村, 高山著 『C によるスプライン関数』 (東京電機大学出版局)[2] C-C. Chang, H-p. Chang, and C-S. Hsue, Chem. Phys. Lett. 217, 486 (1994).
[3] L. Laaksonen, P. Pyykkö, and D. Sundholm, Intern. J. Quantum Chem. 23, 309 (1983).

1Pp041
ネットワーク系における量子アルゴリズム(金沢大理) ○坂本功、山口智也、長尾秀実、西川清【序】量子コンピュータとは、量子力学の原理に基づく全く新しい計算機である。量子コンピュータを実現できれば、従来のコンピュータでは計算時間がかかりすぎて解けない問題
も、高速に解くことができると考えられている。従来のコンピュータでは、ビットは 0 ある
いは 1 のどちらかしか表現できなかったが、量子コンピュータのビット(量子ビット)は 0 と1 の重ね合わせ状態を実現でき、それによって超並列計算を可能にする。N 量子ビットを用いたなら 2 N 通りの状態を実現できる。例えば、2 量子ビットなら 00 , 01 , 10 , 11 の４通りの状態を同時に取り扱うことができる。この重ね合わせ状態と１量子ビットと２量子ビット間の任意のユニタリー変換が実現できれば、量子コンピュータが実現できるといわれている。
１９９４年、Shor による量子コンピュータによる素因数分解アルゴリズムの発見によって、量子コンピュータが全世界で注目を浴びるようになった[1]。以来、ハードウェア実現のための実験的な研究も活発化するようになった。実際には量子コンピュータを実現するために候補となるハードウェアが数多く考えられている。その中で MIT の Issac chuang らによる NMR
量子コンピュータ[2]における７量子ビットの実装が最先端であるが、実用的な数十量子ビットの実現の道筋は現段階では分かっていない。その原因として量子ビットの拡張性が困難であるということが挙げられる。実際に量子コンピュータを実用化するには多くの量子ビットを集積化する必要がある。そこで、現在注目されているのが JST、NTT や NEC によって研究されている電子スピン[3]、量子ドット[4]や超伝導クーパー対[5]などの固体ナノデバイスを用いた量子コンピュータである。固体ナノデバイスでは量子ビットの拡張は易しいが、拡張に伴って量子ゲート数も増えていくのでコヒーレンス時間の短縮や２量子ビット間のユニタリー変換の実現が問題となってくる。その問題を解決する方法として量子誤り訂正や量子ゲートの最適化などが考えられている。そこで、本研究では遺伝的アルゴリズムを用いて固体ナノデバイスにおける量子ゲートの最適化を試みる。
【方法】GA とは系に存在する複数のパラメータを生物界の遺伝の方則を利用して変化させ、ある問題の最適解を求めようという手法である。GA の手順として、コーティング、初期集団の生成し、評価、選択、交叉・突然変異を繰り返し最適解を求める。本研究の場合、
まず固定長の遺伝子の列とアダマール変換や C-NOT などの量子ゲートとを対応づけるコー
ディング作業を行う。次に初期集団を生成するために、乱数を利用して遺伝子配列からなる個体群を生成させ遺伝子列で表される量子回路を計算する。その計算結果が求めたい解か、あるいはその解に近いか、また量子ゲートの数などで個体群を評価の高い順に並べ評価の低い個体群を淘汰する。そして、評価の高い個体群同士の遺伝子を交叉、突然変異させ新しい固体群を先程淘汰した分だけ生成し、また評価するという作業を繰り返し最適解を求める手
法をとる。

１次元型

円型

２次元型

図１

固体ナノデバイスでは、相互作用は隣り合った量子ビット間でしか働かない。よって、遠く離れたところとの相互作用、つまり C ―NOT を実現するには Swap ゲートを用いてその２量子ビット間のユニタリー変換を実現する。しかし多量子ビットにおける C ―NOT を考えた
とき、Swap ゲートの増加により計算ステップ数も増え、デコヒーレンスが効いて量子計算が不可能になる。
これを防ぎ、できるだけ高い確率で解が求められるような最適な量子ゲートを GA を用いて生成する。また、量子ビットの配列についても考える。現在考えられている配列として、
図１のように量子ビットを１列に並べる１次元型が通常用いられるが[6]、本研究では図２のような縦横に並べる２次元型、円型など様々な配置についても考える。これらのモデルで多量子ビットの相互作用を考える時、１次元型では隣り合う量子ビットが２個だが、２次元型
では４個になるので Swap ゲートが少なくてすむことで１次元型より簡単な量子ゲートが生成されると考えられる。また様々な量子ビットの配置における量子ゲートを比較し、ナノデバイスでの量子ビット間の相互作用と量子計算との相関などを検証する。詳細は当日報告する。

参考文献[1] P.W. Shor, SIAM J. Comput,(1997) 1484.
[2] L.M.K. Vandersypen, M. Steffen, G. Breyta, C.S. Yannoni,
M.H sherwood and I.L. Chuang, Nature,(2001) 883.
[3] J.M.Kikkawa, I.P. Smorchkova, N. Samarth and
D.D. Awschalom, Science,(1997) 1284.
[4] H. Kamada, H. Gotoh, J. Temmyo, T. Takagahara and
H. Ando, Phys. Rev. Lett,(2001) 246401.
[5] Y. Nakamura, Yu.A. Pashkin and J.S. Tsai, Nature,(1999) 786.
[6] ALexandre Blais, Phys. Rev. A,(2001)
022312.

1Pp042

第 3 遷移金属 (Hf-Hg) に対する相対論を考慮した電子相関用基底関数(北大院理) ○野呂武司、(青森大) 長内有、(室工大工) 古賀俊勝、(九大院総理工) 三好永作、(苫駒大) 関谷雅弘

【はじめに】原子番号 72 番の Hf から 80 番の Hg までの第 3 遷移金属を含む分子計算
では、相対論の効果とともに電子相関効果が重要であるが、電子相関用の基底関数はほとんど用意されていない。そこで、本研究ではこれらの原子に対してコンパクトでありながら高い記述性を有する電子相関用基底関数を開発した。
【方法】これまで、われわれは原子の高精度の自然軌道 (NO) を再現するように縮約 GTF(CGTF) の係数と軌道指数を決め、高精度の電子相関用基底関数を開発してきた。本研
究でも従来の方法を踏襲し、軌道指数と係数を各軌道量子数 (l) ごとに以下に定義される
汎関数 ∆ を最小にするように決定した。

∆=

Nf
X

Z

nk

d~r|λk (~r) − ψk (~r)|2 w(~r),

(1)

k=1

ここで、 λk (~r) は正確な NOs、nk は占有数、Nf はその個数、w(~r) は重み関数 (1/r2 ) である。また、ψk (~r) =MX

Ci χi (~r),

(2)

i=1

であり、{χi } は CGTF 基底関数、M はその大きさである。χi として、従来は占有軌道すべてを考慮していた。これは、基底関数が全電子計算で利用されることを前提としていたためである。しかし、今回はモデルポテンシャル計算での利用も想定したため、6s、5p、5d、4f の存在だけを考慮した。この結果、得られる基底関数には、今回考慮しなかった
内殻軌道との直交性を記述する節 (ノード) が基底関数に含まれることになる。
NO は藤永等の well-tempered set[1] 拡張した (28s24p20d17f 17g) を用い、各原子の 5dn 6s2 の電子配置での最低状態について原子価電子に関する 1-2 電子励起 CI 計算より求めた。原子積分の計算では、中嶋-平尾 [2] による 3 次の Douglas-Kroll 変換を 1 電子積分について用い相対論の効果を考慮した。基底関数の大きさとしては、[1p1d1f ] と[1s2p2d2f 1g] を、縮約の項数としては、s に対して (5)、p に対して (4)、(3/2)、d に対して (4)、(3/2)、f に対して (5)、(4/2)、g に対して (2) を用いた。全般的に内殻軌道との直交性を保障するために、節のない基底関数に比べて項数が大きくなっている。
【基底関数】表-1 に、開発した基底関数の例として、Os の 1p1d1f を掲げる。
得られた基底関数には、予想通り節の存在が確認できる。モデルポテンシャル MCP ある
いは AIMP では、モデルハミルトニアンの中に内殻軌道の成分をエネルギー的に上方にシフトする演算子があるため、加える基底関数の内殻軌道との直交性を保障する必要がある。したがって、場合によっては直交性を得るためだけに、本来必要のない内殻の軌道を計算の中にあらわに考慮する必要が生じる。今回の基底関数はこの直交性がすでに保障さ表 1: Os の電子相関用縮約ガウス型基底関数

Basis set
[1p1d1f ]lpdf

軌道指数と縮約係数
6.552605 1.182475-.014193.176423
1.475944.355290-.114115.383672
91.724911 24.916424-.016984-.093715.196805-.743108.175543.544315
7.925204-.166347.079337-.359303.076521.205111.841514 .285924.558651 .563594

れているので、このようなことに気を使う必要は全くない。もちろん、他の ECP および
全電子計算でも、問題なく使用できる。
【テスト計算】各原子について、開発した基底関数を用いてテスト計算を行った。占有軌
道は (28s24p20d17f 17g) で求めた HF 軌道を用い、相関軌道として今回開発した [1p1d1f ]、[1s2p2d2f 1g] を用い原子価電子に関する 1-2 電子励起 CI 計算を行った。 比較のために同じ大きさの NO を用い同様の計算も行った。表 2 に、Os の結果を掲げた。
表 2: Os ( 5d6 6s2 (5 D) ) のテスト計算の結果

EHF (a.u.)
ECI (a.u.)
Ecorr (a.u.)-17251.069454-17251.256114-0.186659
NO[1s2p2d2f 1g]-0.178718[1p1d1f ]-0.137974

cGTF (%)-0.177302 ( 99.20 % )-0.137411 ( 99.59 % )今回開発した [1p1d1f ] および [1s2p2d2f 1g] とも、正確な NO と比較して、99% を超える
電子相関エネルギーを与え高い能力を示している。他の原子についても同様の結果を得ている。今回開発した基底関数は、h は含まれないが、[1s2p2d2f 1g] と (28s24p20d17f 17g)との差がすでに 0.007941 a.u. と非常に小さいため、h を含む大きなセットの開発は行わなかった。発表では、モデルポテンシャル計算でのテスト計算および分子系への応用についても述べる予定である。

参考文献[1] (a) S. Huzinaga and B. Miguel, Chem. Phys. Lett. 175, 289 (1990); (b) S. Huzinaga
and M. Klobukowski, Chem. Phys. Lett. 212, 260 (1993).
[2] T. Nakajima and K. Hirao, J. Chem. Phys. 113, 7786 (2002).

1Pp043

誤差関数の計算表式と求値誤差（北教大釧路）小原繁

平成 15 年 7 月 23 日

1

これらの大小関係などにより Fm (T ) と Gm (T ) の値が以下の範囲になることが判かる。

【序】
近年のコンピュータの発達に伴ない分子軌道計算が大型分子系へ適用できる様になってきて、比較的小型の酵素について分子軌道計算が行われ報告される様に

1,
2m + 1q(2m − 1)!! π
0 < Gm (T ) < Sm (T ) =,2(2T )m
T
0 < Fm (T ) < F0 (T ) =なってきた。大型分子系の分子軌道計算を実行するには数値計算の精度についての吟味が不可欠で、特に電子反発積分はその個数が膨大になるにも関らず個々の積分値の有効数字が１２桁程度になっているため計算精度の不足を生じる原因になりやすい。電子反発積分の計算精度を決定している誤差関数を従来よりも高速かつ高精度に求値することを目指し、誤差関数の求値
誤差と計算表式を再検討した。この結果について報告

0 < G0 (T ) <

e−T
.
2T

(8)
(9)
(10)

また、Fm (T ) と Gm (T ) は m や T が充分大きな値のときは以下の値になる。
Fm (T ) → 0, (m → ∞)q
(2m − 1)!! π
Fm (T ) →,2(2T )mT
Gm (T ) → 0.
(T → ∞)(11)(T → ∞)(12)(13)

Fm (T ) と Gm (T ) を T に関して微分することにより次
次節には誤差関数に関する種々の関係式を列挙する。 の関係を得ることができる。

する。

dFm (T ) = −Fm+1 (T ),dTdGm (T ) = −Gm+1 (T ).
dT

第三節には求値のための計算式とそのときの相対誤差、
および、許容できる誤差との関係を導出したものをまとめて記す。

2

Fm (T ) =この節には誤差関数と補誤差関数の定義とこれらの関数が満足する等式を列挙する。
誤差関数 Fm (T )、補誤差関数 Gm (T )、および、和関数 Sm (T ) を下記の様に定義する。

e+ 2T Fm+1 (T ) ,2m + 1ª
1 © −T
Gm (T ) =e
+ (2m − 1)Gm−1 (T ) .
2T

(16)

(17)

この漸化関係を何度か使用すると次式が得られる。

(1)

e−T X
Fm (T ) =am,n + am,N Fm+N (T ),2T(2)

Gm (T ) =

N

Z 1

2

t2m e−T t dt,

Z0 ∞

(18)

n=1

e−T X
bm,n + bm,N +1 Gm−N −1 (T ).
2TN2t2m e−T t dt,

Gm (T ) ≡(15)

部分積分を行なうことにより Fm (T ) と Gm (T ) が次の漸化関係を満足していることが判かる。
© −T
ª
1

【誤差関数に関する等式】

Fm (T ) ≡(14)

1

Sm (T ) ≡ Fm (T ) + Gm (T )q
(2m − 1)!! π=.
2(2T )m
T(T > 0, m = 0, 1, 2, ...)

(3)

被積分関数を吟味することにより Fm (T ) と Gm (T ) に
m や T に関する以下の様な大小関係のあることが判かる。
FM (T ) < Fm (T ),(m < M )

(19)

n=0

(4)

Gm (T ) < GM (T ),(m < M )(5)

Fm (T ) < Fm (t),(t < T )(6)

Gm (T ) < Gm (t).

(t < T )(7)

ここで、am,n と bm,n は
am,0 ≡ 1,
2T
am,n ≡· am,n−1 ,2m + 2n − 1
2T
2T
2T
·
···
,
=
2m + 1 2m + 3
2m + 2n − 1
bm,0 ≡ 1,
2m − 2n + 1
bm,n ≡· bm,n−1 ,2T
2m − 1 2m − 3
2m − 2n + 1
=
·
···
,
2T
2T
2T
(n = 1, 2, . . .)

(20)

(21)
(22)

(23)

である。Fm (T ) を T0 における n 項の Taylor 展開で表現したときの表式は
Fm (T ) =n−1
X(T0 − T )i

i!

i=0

Fm+i (T0 ) + Rn

(24)になる。この式 [(34) 式] から 44 以上の T 値において(33) 式で Fm (T ) を計算すると相対誤差が 10−20 以下になることが判かった。
44 以下の T 値における Fm (T ) の計算には (18) 式の第一項を使うことができる。

(T0 − T )n
Rn =Fm+n (ξ),
n!
ξ ≡ T0 + θ(T − T0 ),(0 < θ < 1)

Fm (T ) ≈

(26)

n
X

Fm (Ti )Ni (T ) + Qn (T ),

(27)

i=1

ここで、Ni (T ) はnY

Ni (T ) =j=1(j6=i)

T − Tj.
Ti − Tj

この場合の相対誤差 ² と許容誤差 σ の関係は²=

am,N Fm+N (T )≤σFm (T )

πn (T ) =(28)
am,N ≤ 2σT eT(−1)n Fm+n (ξ)
,
n!

n
Y

(T − Tj ),(29)

am,n .

(37)

Fm (T ) を T0 における n 項の Taylor 展開で表現した
とき
n−1
X(T0 − T )i
Fm (T ) ≈

(30)

j=1

(31)

を表わす。もし、分点 Ti として、Chebyshev 多項式
Tn (cos θ) = cos nθ の零点µ¶
Tn − T1
T1 + Tn+cos
2
2 cos(π/2n)

N
X
n=1

T1 < ξ < Tn

Ti =

(36)

になる。
この関係式は、予め計算してある誤差関数の値から相対誤差を算出するためには役立つが、許容誤差以内
で誤差関数の値を計算するには N をどのような値にすればよいかを決めるためには役立ちにくい。このための関係式を相対誤差 ² の表式 [(36) 式] に (8) 式と (18)式を代入して変形することにより以下のように導出することができた。

を表わし、Qn (T ) は剰余項
Qn (T ) = πn (T )

(35)

n=1

(25)

を表わす。また、n 個の T の値 (T1 , T2 , . . . , Ti , . . . , Tn )における誤差関数の値 Fm (Ti ) を使って Fm (T ) の Lagrange 補間の表式を得ることができる。
Fm (T ) =

e−T X
am,n
2T
N

になる。ここで、Rn は剰余項でありn − k + 1/2
πnFm+i (T0 )

(38)

|Rn |
|T0 − T |n Fm+n (ξ)=
Fm (T )
n!Fm (T )

(39)

i!

i=0

相対誤差 ² は²=になる。Taylor 展開における T0 のグリット幅を w とし
T = T0 + w, (32)(40)

相対誤差 ² が許容誤差 σ 以内であるという条件を満たすグリット幅 w の表式として次式を得ることができた。

k = 1, 2, . . . , n

を 選 ん だ 場 合 、こ の Lagrange 補 間 の こ と を 特 に|w|n ≤

Chebyshev 補間と呼ぶ。

σn!Fm (T0 )Fm+n (T0 − |w|)

(41)

Lagrange 補間式で Fm (T ) を計算したときの3

Fm (T ) ≈【計算表式と相対誤差】

nX

Fm (Ti )Ni (T )

(42)

i=1

誤差関数に解析的な計算表式が無いので誤差関数の相対誤差 ² は値を計算するには何らかの近似式を使うことになる。

²=

前節の等式に基づいた計算表式とその相対誤差 ²、お|Qn ||πn (T )|Fm+n (ξ)=
Fm (T )
n!Fm (T )(43)

よび、許容できる誤差 σ との関係をこの節に列挙する。 になる。π (T ) の大まかな見積りとして区間幅 (T −T )nn1(19) 式の第一項により Gm (T ) を計算し、この値を の n 乗を用いると
Sm (T ) から差し引くことにより Fm (T ) を計算するこ|Tn − T1 |n Fm+n (ξ)とができる。
(44)
²<
e−T X
bm,n
Fm (T ) ≈ Sm (T ) −2T

n!Fm (T )

m−1

(33)

この場合の相対誤差 ² と許容誤差 σ の関係は²=

(2m − 1)!!e−T
bm,m G0 (T )<≤ σ.
Fm (T )
(2T )m+1 Fm (T )右辺は (39) 式の右辺と同形になる。区間幅 (Tn − T1 )を w で表わすと、許容誤差と w の関係は Taylor 展開

n=0の時と同じ関係になることが判かった。
以上の計算表式と相対誤差の詳細については当日発
(34)

表する。

1Pp044
スペクトル分解法による Schrödinger 方程式の数値解法(II)(名工大工 1、名工大院･おもひ領域 2) ○滝本理人 1、志田典弘 2【序】振動や回転等の原子核の運動に対して量子学的に固有解を求める方法は、Numerov 法、有限要素法、DVR 法等数多くの方法が知られているが、一般的かつ複雑な問題にまで対応できるような方法は現在まで確立されていない。そこで我々は、複雑な系や高エネルギー
状態等の多様な問題にも対応可能な方法としてスペクトル分解法に注目し、この方法を洗練し実用化のレベルにまで高めるための工夫を行ってきた。特に昨年は、高エネルギー状態を精度良く求めることに焦点を当て、初期波動関数の選択やノイズの除去等に関する工夫を行い、その成果を本討論会において報告した。今回の発表では、多次元への拡張とそれに伴う計算時間、計算量の軽減に対する工夫について報告する。
【計算方法】
スペクトル分解法は、時間依存の波動関数が、定常状態の固有関数と、固有値に対応した位相因子との積の重ね合わせであることに注目し、フーリエ変換で位相因子の違いを基に分解することによって固有値と固有関数を求める数値解法である。その具体的な計算方法の概略は次のようである。まず波動関数を Split Operator を用いて時間発展させることで波動関数の時系列を得る。
i t
i t
Ψ (x, y, z , t 0 +Δ t ) = exp ∆ ∇ 2 × exp(−i ∆ tV )× exp ∆ ∇ 2 × Ψ (x, y, z , t 0 )× O[(∆t ) 3 ] 4M 
 4M 

①

ここで Split Operator は波動関数の微小時間変化における時間発展を、右辺第一，三項によって表される自由粒子としての伝播と第二項によって表されるポテンシャルの寄与とに分離して行うものであり、これは波動関数の表示形式をそれぞれ運動量表示と座標表示にフーリエ変換で変換することで可能となる。得られた波動関数の自己相関関数を以下の式により求める。

ρ (t ) = Ψ (x, y , z ,0 ) Ψ (x, y , z , t )②

そしてこの自己相関関数に対し、Wiener-khinchin の定理に従ってフーリエ変換を行いパワースペクトルを求める。
+∞

1
S (ω ) =ρ (t ) exp(− iω ⋅ t )dt
2π −∫∞③

なお、実際に③式を評価するときには Window Function を付加し、計算精度を安定させる。

スペクトル分解法では、ここで得られたパワースペクトルのピークの位置より固有値を求め、更にこの固有値に対応した位相因子を時間依存の波動関数に掛けて時間で積分することにより固有関数を求める。
Tf

Ψn ( x, y , z ) = ∫ dt ∑ C i Φ i ( x, y , z ) exp[i (E n − E i ) ⋅ t ]0④i

この方法の大きな利点として、全ての固有状態を同等に扱うために理論上高エネルギー状態においても精度の低下がないこと、多次元の問題への拡張が比較的容易であること
等が挙げられる。しかしながら、この方法を実際に化学の問題に応用するためには、計算精度や計算時間といった数値計算上のいくつかの問題を解決する必要がある。このうち高エネルギー状態を求める際に発生する計算精度の問題については、序で述べたように改善法を考案した。そこで、今回の講演では多次元とそれに関連する問題点に対する改良について報告する。スペクトル分解法においては多次元化の際に次元数の増加とともに計算量が飛躍的に増大し、その増加量はおおよそ次元数のべき乗となる。従ってこのままでは理論上何の制約がなくとも２～３次元の計算で破綻がきてしまう。そこで我々は Split
Operator による波動関数の時間発展をすべて運動量空間で行うことにより計算量の軽減を図る方法を考案した。又、この方法は計算時間だけでなく必要とされる計算機リソース（メモリやディスク）も次元数のべき乗で増大する。そこで我々はこの点についても計算アルゴリズムを洗練することにより解決した。
【計算結果】
多次元の計算のプロトタイプとして２次元のモデル系の計算結果を示す。これはHamiltonian が(1)式で与えられる、
プロトン移動反応に類似したモデル系である。表１は、得られたモデル系のエネルギーである。表１より、このような２次元の計算でも系のエネ
ルギーを正確に再現できることがわかる。この例からも、スペクトル分解法が分子振動等を計算するための,多次元の Schrödinger 方程式の数値解法として有用な方法であることがわかる。又、現実系に即したより複雑な問題への応用は、方法論についての詳細と合わせ当日報告の予定で
ある

Ĥ =

h2 ∂2 h2 ∂2++ V ( x, y )2m ∂x 2 2m ∂y 2

(1)

V ( x ) = 0.01x 4 − 0.02 x 2 + 0.01

+ 0.0005 y 2

m = 1000 (a.u)

表１， モデル系の計算結果

nE n 計算値 × 10

3

(a.u)

E n 理論値 × 103

１4.6136

4.6136２

4.73874.7387３

5.6136

5.6136４

5.73875.7387５6.6136

6.6136６

6.73876.7387

1Pp045

局所制御理論による振電励起制御レーザー場の設計
慶大院・理工菅原道彦【１】序
近年のレーザー技術によってもたらされる高強度フェムト・アト秒領域の超短パルスは、分子振動の時間スケール領域に達しており、ピーク強度は分子を強制振動させるほどの出力を併せ持っている。これらの事実は、分子と高強度・非定常レーザー場とのコヒーレントな相互作用を利用して望みの量子力学的非平衡状態を作り出し、積極的に分子量子状態を目的の状態に遷移させるといった「量子ダイナミクスの制御」の可能性を示唆している。すなわち、レーザー場を精密に設計することにより電子移動、イオン化、化学反応などの様々な物理現象に対するミクロスコピックな観点からの量子制御が期待される。我々は量子系の運動方程式の形を考慮し各時刻での系の情報のみを用いて（時間領域で）局所的に制御場を設計する
局所制御理論に注目し研究を行ってきた。この理論の大きな特徴として、最適制御理論に基づいた理論において不可避である数値的な扱いが困難な二点境界値問題を回避できるということが挙げられる。これまで、我々はこの一般的局所制御理論を準位分布の制御、波束動
力学の制御（終状態指定、トラッキング問題等）に適用してきた。本研究では、局所制御理論の適用範囲を紫外レーザー領域にも広げ、電子励起過程を利用した振動波束制御を可能とするレーザー場の設計を目指す。この際、フロケ描像の下で系を定式化することによりパルスのエンベロープのみを最適化の対象とする。これにより、電子遷移に対応するキャリア振動数を一定に保った形で制御場の設計ができるため、実際のレーザー生成においてより現実的であるレーザー場を提示できる。
【２】一般理論
高強度・非定常レーザー場による２つの電子状態間の振電励起を考慮するためにフロケ描像を採用する。回転波近似の下で系のシュレディンガー方程式はd ψ g (t ) i  H + hω = −  gdt ψ e (t ) h  µ ⋅ u (t )()µ ⋅ u (t )  ψ g (t ) ⋅H e  ψ e (t ) 

,

(1)

i
≡ − Hˆ 0 + Vˆu (t ) ⋅ Ψ (t )hと書くことができる。但し、 H g 、 H e 、 µ はそれぞれ電子基底状態、励起状態のハミルト
ニアン及び遷移モーメント演算子である。また、ω はパルスエンベロープ u(t)のレーザーパ
ルスの中心角振動数を示している。
ここで、現在の状態が目的の量子状態とどの程度合致しているかを評価する量として()

y (t ) = y Oˆ 1 (t ) , Oˆ 2 (t ) ,L , Oˆ N (t ) を導入する。 Oˆ i (t )(i = 1,2,L, N ) は i 番目の観測量に対する演算子 Oˆ i (t ) の期待値であり、演算子 Ôi (t ) には運動方程式[]dOˆ i (t )= i Oˆ i (t ), Hˆ 0 ,dt(2)を要請する。時刻 t に対して y (t ) を微分し、 シュレディンガー方程式(1)式及び(2)式を考慮することにより、評価指数の時間微分はN∂y (t )dy= 2∑
Im Oˆ i (t ) Vˆ u (t ) ,ˆ
dt
i =1 ∂ Oi (t )と表される。ここで、制御レーザー場エンベロープを(3)

∂y (t )Im Oˆ i (t ) Vˆ ,
ˆi =1 ∂ O (t )Nu (t ) = A0 ∑(4)i

とおくことにより、 y (t ) の単調増加を保証するレーザー場を設計することが出来る。
【３】応用
モデル系として電子基底状態に調和振動子型、励起状態としてモース型の断熱ポ
テンシャルを持つ分子を考える。電子基底する振動数のレーザーを照射した場合フ
ロケ描像における二つの断熱ポテンシャ
ルは図１に示すように交差する。今回は制e

8

potential energy

g 、励起状態間 e のエネルギー差に相当

10

御量演算子 Oˆ i (t ) として電子励起状態上の

6

4

2

振電状態への射影演算子を採用すること
により、電子励起を促進するパルスの設計

g-5を行う。制御量演算子 Oˆ i (t ) は系のハミル

5

0

図１

10

15

20

25

x

フロケ描像下の断熱ポテンシャル

トニアンと交換可能であるため運動方程式(2)式の右辺は０となり、結果として Oˆ i (t ) は時間
発展しない。このため、レーザー場は(4)式に Oˆ i (t ) の初期状態演算子をそのまま代入することによって得られる。シュレディンガー方程式の時間発展は４次のルンゲ・クッタ法で数値的に解いている。
図２に設計されたパルスエンベロープの下での基底状態、励起状態の準位分布の時間変化を示す。設計されたレーザー場の下で、準位分布が電子基底状態から励起状態へほぼ１００％移動している様子が示されている。図３は得られたレーザーパルスのエンベロープの時間変化を示している。得られたエンベロープによると、制御レーザー場が二つのパルス列として設計されていることがわかる。図２には二つのパルスが照射される時間帯に対応し、分布が二段階の過程を経て励起されている様子が示されている。
振動数が時間変化するようなチャ－プパルスのエンベロープ設計も本理論で取り扱うことが出来る。また、二つの独立な系をハミルトニアンに組み込んで制御場を設計することにより、同位体分子の選択的励起する制御場を設計することも可能である。これらの結果については当日発表する予定である。
1

g

e

1.5

amplitude

population

0.8

0.6

0.4

1

0.5

0.2

0
20

40

60

80

100

0

time

図２

準位分布の時間変化【参考文献】
M. Sugawara, Chem. Phys. Lett., 358, 290 (2002).
M. Sugawara, J. Chem. Phys., 118, 6784 (2003).

20

40

60

80

time

図３

設計されたエンベロープの時間変化

100

1Pp046
フラグメント MO 法を用いた量子・量子および量子・古典融合法（産総研・計算科学）○北浦和夫、Dmitri G. Fedorov、古明地勇人、上林正巳【序】我々は、数千原子からなる分子・分子集合体の ab initio 計算を可能とする近似的計算方法（フラグメント MO 法；FMO 法）を開発している。1) FMO 法は、分子・分子集合体を小さなフラグメ
ントに分割し、フラグメント（モノマー）とフラグメントペア（ダイマー）について MO 計算を行うだけで、全系の全エネルギーとプロパティが計算できる方法である。この方法は、通常の ab initio 計
算の結果を非常に精度よく再現する。また、高効率並列計算が可能なことも一つの特徴である。現在開発中の FMO 法を組み込んだ Gamess によると、128 台の PC で並列化効率 80〜90%を達成している。
このシステム（CPU は Pentium3 の 1GHz、 FastEthernet 接続）による lysozyme（約 2000 原子系）の 1
点計算（FMO-HF/STO-3G）の elapsed time は約 30 分である。
FMO 法を用いると数千から数万原子系の HF/double-ζ程度の 1 点計算は容易に行える段階にあるが、実用的には、より大きな基底関数で電子相関まで考慮して、さらに構造最適化や分子動力学シミュレーションを行うことを考えると、より一層高速な計算方法が必要となる。数千原子系全体を abinitio 法で扱って、このような計算が行える目途は今のところたっていない。現実的なアプローチとしては、ONIOM 法 2)のような量子・量子、量子・古典ハイブリッド法が有用である。系の重要な部
分は電子相関レベルで、その他の部分は HF レベルで、さらには大部分のあまり重要でない部分は古典力場で扱うことにより、大規模系の高速な計算が可能となる。本発表では、分子を分割して計算するという FMO 法の特徴を生かした ONIOM 法様のハイブリッド法について報告する。
【方法】まず、FMO 法の概要について述べる。FMO 法は、分子や分子集合体をいくつかの小さなフラグメントに分割し、フラグメント（モノマーと呼ぶ）とフラグメントペア（ダイマーと呼ぶ）について ab initio MO 計算を行うだけで、全系のエネルギーとプロパティを計算する方法である。一般的にフラグメントサイズを小さくするとエラー（通常の ab initio MO 計算結果との差）が大きくなり、
10~20 原子程度（ポリペプチドでは、2 アミノ酸残基程度）以上にとると、全エネルギーのエラーは数 kcal/mol 程度となる。モノマーとダイマーの計算は、通常の 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)

i

x はモノマー(x=I)またはダイマー(x=IJ)を示す。一電子ハミルトニアンに加えた射影演算子項は、フラグメント境界にある原子の基底関数をフラグメントに割り振る働きをする。また、V は静電ポテンシ
ャルである。通常の HF 法と異なるのは、例えば、モノマー I の計算には、I 以外のすべてのモノマーからの静電ポテンシャルが含まれることである。同様に、ダイマー IJ の計算には、IJ 以外の全てのモノマーからの静電ポテンシャルが含まれる。分子の全エネルギーは、モノマーとダイマーの全エネル
ギー、EI と EIJ を用いて、次式により計算する。

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

(3)

I

Exに含まれる静電相互作用エネルギーをまとめると、式(3)は次ぎのようにフラグメント間相互作用エ~
ネルギー, ∆E IJ ,を用いて書きかえることができる。3)~~E = ∑ E I + ∑ ∆E IJ
I

(4)

I >J

FMO 法を用いた融合法では、FMO 法のフラグメント境界に合わせて領域を設定することにより、波動関数や基底関数を混在させることが容易にできる。また、最も低

いる。

 ~
~ ( L1) E =  ∑ EIFMO( L1) + ∑ ∑ ∆EIIFMO
'

I
L
1
I
L
1
I
'
(
I
)
L
1
∈
∈
>
∈


 ~
~ ( L2) ~ + ∑ ∑ ∆EJIFMO( L2)+  ∑ EJFMO( L2) + ∑ ∑ ∆EJJFMO'J∈L2 J '( > J )∈L2
 J∈L2
 J∈L2 I∈L1

MM 
MM
MM

+  ∑ EKMM + ∑ ∑ ∆EKK
'  + ∑ ∑ ∆EKJ + ∑ ∑ ∆EKI
K∈L3 K '( >K )∈L3
K∈L3 I∈L1
 K∈L3
 K∈L3 J∈L2

fI’

fJ’

レベルの領域に分子力場を採用して、QM/MM 法として用いることができる。全系のエネルギー計算には次式を用

Int1-1

Level1fK’

Int2-2 fJ

Int3-3 f
K

fI
Level2Int2-1Int3-2Int3-1

Level3

図 1 領域分割と領域間の相互作
用エネルギー計算の概念図(5)

現在、汎用なコードを Gamess に組込んでいるところであるが、現バージョンの量子・量子融合法
では、基底関数は任意のものを混在させることができ、波動関数は、HF をはじめ、DFT、MP2、MCSCFを混在させることができる。さらに、IMOMM4)タイプの QM/MM 計算が行える。QM/MM については、
領域境界を水素原子でキャップしなくてすむ方法を開発中である。
【結果】HF レベルで基底関数を混在させて、lysozyme と tri-[Nacetyl-glucosamine]複合体(図 2)における、タンパク質とリガンドの相互作用エネルギーを計算した。リガンドとリガンドを取り巻く 16
アミノ酸残基を FMO-HF/6-31G で、その他の部分は FMO-HF/STO-3G

HF/STO-３Gを適用した。スーパーモレキュール計算で求めた分子間相互作用エ
ネルギーは−88.9kcal/mol となり、全系を FMO-HF/6-31G で計算して得られた相互作用エネルギー−87.6kcal/mol と極めて近い値が得

HF/6-31Gられた。
（約 1kcal/mol の過大評価。
ちなみに、全系を FMO-HF/STO-3Gで計算した場合の相互作用エネルギーは−53.4kcal/mol）これにより融合法の精度に関する検証はできたが、問題点としては、基底を混合した計算では収束性が悪くなり、モノマー SCF の収束に約 100 サ
イクルかかったため、計算時間の大幅な短縮にはならなかった。
現在、この点について改良を進めている。

図 2 lysozyme と NAG3 の複合体。
（PDB:1HEW）

1) K.Kitaura et al., Chem.Phys.Lett.,312,319-324(1999)
2) M.Svensson, S.Humbel, R.J.Froese, T.Matsubara, S.Sieber, K.Morokuma, J.Phys.Chem., 100(1996)19357
3) T.Nakano et al., Chem.Phys.Lett., 351, 475-480(2002)
4) F.Maseras, K.Morokuma, J.Comput.Chem.,16(1995)1170