4P101
シクロデキストリンと 胆汁酸・胆汁塩の複合体についての理論的研究（お茶大院人間文化創成科学）姚 嵐、森 幸恵、鷹野 景子

Theoretical study on complexes of cyclodextrins with bile acids
and bile salts(Ochanomizu Univ.) Lan Yao, Yukie Mori, Keiko Takano
Introduction
Cyclodextrins (-CD, -CD, and-CD) are cyclic oligosaccharides consisting of six,seven and eight -(1,4)-linked D-glucose units, and they are good host molecules toaccommodate hydrophobic guest molecules inside their cavity. According to the previous
calculation in the gas phase, the lowest energy conformation is the closed conformation (Fig.
1a), while another energy minimum is the open conformation, which resembles the structures
observed in the solid state (Fig. 1b).1 It is reported that CDs form 1:1 inclusion complexes
with bile salts such as cholate (CA-) and deoxycholate (DCA-) in aqueous solution. In thisstudy, the intermolecular interactions have been investigated for 1:1 complexes of CDs withcholic acid (CA, Fig. 2a), deoxycholic acid (DCA, Fig. 2b) and their anions (CA-, DCA-) inthe gas phase and in aqueous solution by means of the DFT calculations.

(a)

(b)

(a) CA(b) DCA

Fig. 1. Top and side views of -CD in the open

Fig. 2. The structural formulas of cholic acid(a) and closed (b) conformations.

(CA) and deoxycholic acid (DCA)

Computational methods
The geometries of the CD/CA inclusion complexes were optimized by the DFT-D(B97-D) and DFT (M06-2X) methods with the standard 6-31G(d) basis set. The strengths ofinteractions were quantified by the basis set superposition error (BSSE)-corrected binding
energies (Ebinding) both in the gas phase and in aqueous solution. The polarizable continuum
model (PCM) was applied for the aqueous solution. The strain energy (Estrain) was also
evaluated. The intermolecular hydrogen bonds and van der Waals interactions between CA
and CD were investigated by inspection of the optimized geometries

Results and Discussion
According to the results of the calculations by the B97-D and M06-2X methods, theopen/top1 and the open/mid2 are energetically favored configurations for the -CD/CA
complex in the gas phase, while open/bottom1 and open/top1 are favorable for the -CD/CA
complex. The structures are schematically shown in Fig.3, in which the arrow indicates the
direction from the carboxylic group toward the A-ring of CA (defined as direction 1,
otherwise is direction 2). The diameter of-CD (dup =7.13Å, ddown = 6.71 Å in the open
conformation) prohibits the steroid skeleton of the bile acid or bile salt (width = 6.20 Å forCA) from entering into the cavity of -CD, and only the side chain can enter into it. Suchgeometrical features are in good agreement with those suggested by the experimental study.2

Fig. 3. Schematic drawings of the stable configurations of the CD/CA complexes.
(Top, mid, and bottom are defined by the complexed position of the main part of guest molecule.)

Three hydrogen bonds are seen between -CD and CA, suggesting that the H-bonds andvan der Waals interactions play a decisive role in the complexation in the gas phase.
The BSSE-corrected binding enthalpies were calculated by the M06-2X/6-31G(d) method
for the -CD complexes in aqueous solution (Table 1). The present calculations overestimate
the absolute values of Hbinding, but the difference between CA- and DCA- is quite small, which
agrees with the experimental observation. Such small difference between CA- and DCAsuggests that the hydroxyl group at C7 only slightly affects the complexation with -CD.
Table 1. BSSE-corrected binding energies and binding enthalpies /kcal·mol-1 in aqueous solution
Guest

Ebinding

Hbinding

Experimental Hbinding2

CA

-22.5

-19.8

-----

-

-25.0

-22.9

-5.5

DCA

-22.9

-21.1

-----

-

-25.0

-22.5

-6.2

CA

DCA

1. Karpfen et al. Monatsh. Chem. 2008, 139, 363.
2. Liu Y et al. J. Org. Chem. 2007, 72, 8227.

4P102
離散時間マスター方程式を用いた不均一媒体中での拡散方程式の研究
豊橋技術科学大学大学院 情報・知能工学専攻 関野研究室
石井智茂，墨智成，関野秀男
A discrete-time master equation study on diffusion process in inhomogeneous medium
Toyohashi University of Technology Computer Science and Engineering Sekino laboratory
Ishii Tomoshige, Sumi Tomonari, Sekino Hideo
背景1

流動モザイクモデルの提唱により，細胞膜を構成する脂質二重膜は，二次元の流体を形成し，膜タンパク質はモザイク状に埋め込まれる形で，この流体中を動き回りながら，細胞間の情報伝達や，栄養素の摂取および代謝物の排泄等の役割を果たしていると考えられていた．ところが近年，脂質ラフトの提唱により，特定の場所に膜貫通蛋白質や脂質を集合／会合させ，局在化させることで，分子の組織化を行い，様々な機能をもつ膜貫通蛋白質の集合／配列構造を細胞膜上に構築することにより，上記の働きを効率的に行い，さらに外部環境の変化に応じて，細胞膜分子の組織化をきわめて動的に変え，構造体は生成・消滅を繰り返しているのではないかと考えられる様になった．
以上を踏まえ本研究では，この細胞膜の不均一性を空間に依存した拡散係数として考慮することにより，細胞膜の不均一性に起因した膜貫通蛋白質の集合／拡散運動への影響を単純なモデルを用いて解析することを最終的な目標として研究を進める．その準備として，Fick の法則と連続の式から導かれる拡散方程式を出発点にして，不均一媒体中における拡散過程を記述する離散時間マスター方程式を導出し,空間に依存した拡散係数の下でのシミュレーションを行った．

理論

2

本研究ではランダムウォークを解析するために，Fick の法則である式(1)と連続の式(2)から離散時間マスター方程式(3)を導いた．この方程式は解析的に確率分布の保存を満たす．式(1)における𝐷(𝑥, 𝑦)は格子上の座標(𝑥, 𝑦)における拡散粒子の拡散係数を表す．𝜑(𝑥, 𝑦)は格子上の座標(𝑥, 𝑦)での離散時間 t における拡散粒子の確率分布を示す．
𝐽(𝑥, 𝑦, 𝑡) = −𝐷(𝑥, 𝑦)∇𝜑(𝑥, 𝑦, 𝑡)(1)

𝜕𝜑(𝑥,𝑦,𝑡)(2)𝜕𝑡= −∇ ∙ 𝐽(𝑥, 𝑦, 𝑡)𝜑(𝑥, 𝑦, 𝑡 + ∆𝑡) = {1 −−3𝐷(𝑥 + ∆𝑥, 𝑦) + 𝐷(𝑥, 𝑦) 𝐷(𝑥 − ∆𝑥, 𝑦) + 𝐷(𝑥, 𝑦)
−
10𝐷0
10𝐷0

𝐷(𝑥, 𝑦 + ∆𝑦) + 𝐷(𝑥, 𝑦) 𝐷(𝑥, 𝑦 − ∆𝑦) + 𝐷(𝑥, 𝑦)−} 𝜑(𝑥, 𝑦, 𝑡)
10𝐷0
10𝐷0

+

𝐷(𝑥 + ∆𝑥, 𝑦) + 𝐷(𝑥, 𝑦)𝐷(𝑥 − ∆𝑥, 𝑦) + 𝐷(𝑥, 𝑦)𝜑(𝑥 + ∆𝑥, 𝑦, 𝑡) +𝜑(𝑥 − ∆𝑥, 𝑦, 𝑡)
10𝐷0
10𝐷0

+

𝐷(𝑥, 𝑦 + ∆𝑦) + 𝐷(𝑥, 𝑦)𝐷(𝑥, 𝑦 − ∆𝑦) + 𝐷(𝑥, 𝑦)𝜑(𝑥, 𝑦 + ∆𝑦, 𝑡) +𝜑(𝑥, 𝑦 − ∆𝑦, 𝑡)
10𝐷0
10𝐷0

(3)

モデル

本研究では，上記の離散時間マスター方程式を用いてシミュレーションを行い，不均一媒体中での粒子拡散の基本的な性質について解析した．空間に依存した拡散係数𝐷(𝑥, 𝑦)に対して，単純な周期関数の式(4)2𝜋𝑥2𝜋𝑦𝐷(𝑥, 𝑦) = 𝐷0 × {1 − 𝐴𝑥 × sin () − 𝐴𝑦 × sin (
)}
𝐿𝑥
𝐿𝑦

(4)

を仮定したモデルを採用した．ここで Ax， Ay は拡散係数の変化率 Lx， Ly は 1 周期の長さをそれぞれ表す．変化率 Ax， Ay が共にゼロの場合，𝐷(𝑥, 𝑦) = 𝐷0 となることから，均一媒体中での自由拡散に帰着する．すなわち Ax, Ay
および Lx, Ly はそれぞれ不均一性の強さおよび周期を示している．また，𝐷(𝑥, 𝑦)の空間平均値は𝐷0 と一致する事から，本モデルにより，媒体の不均一さに起因した通常拡散係数への影響を調べる事が出来る．

4

結果及び考察𝐷(𝑥, 𝑦)の周期長を 12，24，36 に固定した場合の，ステップ数 400 における平均二乗変位から見積もられる通常拡散係数の変化率依存性を Fig.1 に示す．Fig.1 から分かる様に，𝐷(𝑥, 𝑦)の変化率 Ax, Ay が大きくなるに従い，通常

拡散係数が均一系(Ax,=Ay=0)での自由拡散係数𝐷0 = 0.8より小さくなっていく事が分かる．また，変化率の各値に沿ってみた場合，周期長が長くなるに従い，通常拡散係数が減っていることが分かる．よって Fig.1 から，𝐷(𝑥, 𝑦)の変化率と周期長の両方が，通常拡散係数に影響を与える可能性が示唆される．しかしながらこれらの結果は，ステップ数が少ない事による人為的結果である可能性を含むので，以下でその検証を行う。

ステップ数400
0.805
通
0.8
常 0.795
拡 0.79散 0.785
係 0.78
数
0.775

D

周期長12
周期長24
周期長36

0

0.05

0.1

0.15

0.2

変化率Ax, Ay

Fig.1 ステップ数 400 における各周期長に対する(通常)拡散係数の変化率依存性
ステップ数 400 において通常拡散にいたっているか確認するために，ステップ数 2400 における平均二乗変位から見積もられる通常拡散係数の変化率依存性(周期長 12 の場合)を Fig.2 に示す．

周期長12

0.805

通
常
0.8
拡
散 0.795
係
数
0.79

D

ステップ数400
ステップ数2400
0

0.05

0.1

0.15

0.2

変化率Ax, Ay

Fig.2 周期長 12 における各ステップ数に対する(通常)拡散係数の変化率依存性
Fig.2 から分かるように，ステップ数 400 とステップ数 2400 の結果が完全に一致した事から，周期長 12 の場合、400 ステップでも十分に通常拡散に達していると言える．この事から，拡散係数𝐷(𝑥, 𝑦)の不均一性の強さ(Ax, Ay)が増加するに従って，通常拡散係数は均一系の自由拡散係数より小さくなっていく事が分かる．
次に，通常拡散係数が周期長に依存しているかを長いステップ数のシミュレーションによって確かめるために，
ステップ数 2400 における通常拡散係数の変化率依存性を Fig.3 に示す．

ステップ数2400
0.805
通
常
0.8
拡
散 0.795
係
数 0.79

D

周期長12
周期長24
周期長36

0

0.05

0.1

0.15

0.2

変化率Ax, Ay

Fig.3 ステップ数 2400 における各周期長に対する(通常)拡散係数の変化率依存性
Fig.3 から分かる様に，ステップ数 2400 では周期長の変化による通常拡散係数の変化ほとんど見られず，周期長
24 および 36 の曲線は，十分に定常状態に達している周期長 12 の曲線に重なる傾向が伺える．このことから Fig.1
で示した，周期長 24 および 36 の結果は，400 ステップではステップ数が足りていないため，通常拡散に至っていない段階での結果であると言える．いずれにせよ，Fig.3 から，定常状態に至った際に通常拡散係数に対する𝐷(𝑥, 𝑦)の周期長依存性は極めて小さいか，あるいは全くない可能性が示唆される．発表では，さらに長いシミュレーシ
ョンによる結果との比較を報告する．今回示した結果は全て，x 方向と y 方向の不均一性を等価に扱った場合のみの結果であったが，今後は異方性向を考慮した場合も含め解析を進める予定である．

4P103
フラグメント密度汎関数理論による
β ア ド レ ナ リ ン 受 容 体 -リ ガ ン ド 相 互 作 用 の 解 析（安田女子大・薬）下堂

靖代、柴田

瑩

Interaction Analyses on β adrenergic receptor-ligand complexes
with Fragment Density Functional Theory(Yasuda Women’ s University, Faculty of Pharmacy)
Yasuyo Shimodo, Akira Shibata【序】近 年 ， 人 体 に 存 在 す る G タ ン パ ク 質 共 役 受 容 体 (GPCR)と し て ， 初 め て

β 2 ア ド レ ナ リ ン 受 容 体 の 三 次 元 構 造 が 明 ら か に な っ た [ 1 ] ．G P C R は 薬 物 標 的 タ
ンパク質であることから，この解明により新薬・改良薬開発の効率化が期待される．論理的な薬物設計においては，標的タンパク質の立体構造情報だけでなく，リガンド化合物との相互作用基盤について理解することも不可欠である．
本 研 究 で は ，計 算 科 学 的 な 手 法 を 用 い て ，β ア ド レ ナ リ ン 受 容 体 - リ ガ ン ド 複 合体における結合部位および作用部位の解明を目指している．
【計算方法】

β 受 容 体 に は β 1 ，β 2 ，β 3 の サ ブ タ イ プ が あ る が ，β 1 お よ び

β 3 受 容 体 の 三 次 元 構 造 は 明 ら か で な い た め ，β 2 受 容 体 の 結 晶 構 造 を も と に ホ
モロジーモデリングにより予測した．これらの各サブタイプと非選択的な作動薬（ n o r a d r e n a l i n e , a d r e n a l i n e , i s o p r o t e r e n o l ）お よ び 拮 抗 薬（ p r o p r a n o l o l ,pin, pindolol） と の 複 合 体 を モ デ リ ン グ し ， 310K で 5ns の 分 子 動 力 学 シ ミ ュレーションを行った．これらのシミュレーションでのスナップショットを抽出し，フラグメント密度汎関数理論による電子状態計算から，フラグメント間相互 作 用 エ ネ ル ギ ー Δ EIJ を 求 め た ．
【結果】受容体のアミノ酸残基とリガンドの相互作用エネルギーに相当する

Δ E I J を 解 析 し た と こ ろ ，β 1 ，β 2 ，β 3 受 容 体 の そ れ ぞ れ A S P 1 3 8 ，A S P 1 1 3 ，A S P 1 1 7に お い て ，い ず れ の リ ガ ン ド と も 強 い 引 力 的 な 相 互 作 用 が み ら れ た ．ま た ，β 1 ，
β 2， β 3 受 容 体 の そ れ ぞ れ SER232， SER207， SER212 に お い て は ， 作 動 薬 と 引
力的な相互作用がみられたのに対し，拮抗薬との相互作用はみられなかった．
さ ら に ，こ れ ら の S E R と リ ガ ン ド の Δ E I J に は E C 5 0 と の 対 応 が み ら れ た ．し た がっ て ，β 受 容 体 の 主 な 結 合 部 位 と 作 用 部 位 は ，そ れ ぞ れ 同 位 の A S P と S E R あ り ，
これらは複合体構造上も合理的であることを明らかにした．
[ 1 ] V . C h e r e z o v , et a l . , S c i e n c e 3 1 8 ( 2 0 0 7 ) 1 2 5 8 .

4P-104
経路積分法による周囲環境を考慮した核酸塩基対における水素結合の解析（横市大院・生命ナノ 1, 原子力機構 2） 大道雅史 1, 志賀基之 2, 川島雪生 1, 立川仁典 1

Nuclear quantum effect on the hydrogen-bonded structure of nucleobase pairs
under the surrounding environment with path integral method(Yokohama-city Univ. 1, Japan Atomic Energy Agency2) Masashi Daido1,
Motoyuki Shiga2, Yukio Kawashima1, Masanori Tachikawa1【序論】 核酸塩基対間の水素結合は、DNA の二重らせん構造の維持や遺伝情報の伝達に関わる重要な分子間相互作用
で あ る [1,2] 。 そ の た め 、 多 く の 構 造 異 性 体 と と も に
Watson-Crick 型塩基対の詳細な構造を得るために、実験と理論計算の双方から多くの研究が行われている[3,4]。これまでの塩基対の理論的解析のほとんどは、ab initio 分子軌道法

Guanine-Cytosine (G-C) pairや古典的な分子動力学法によるものであるが、水素結合系のより定量的な記述のためには核自身の量子効果が必須である
ことが知られている[5]。そこで、本研究では温度効果および核の量子効果をともに考慮できる経路積分ハイブリッドモンテカ
ルロ法[6]を用いて、図１に示した Watson-Crick 型の核酸塩基対を計算し[7]、温度変化による塩基対構造における核の量子効果の影響を考察した。また、主成分分析による塩

Adenine-Thymine (A-T) pairFig.1 Structure and atomic labelingof Watson-Crick type base pairs

基対全体の揺らぎも解析した。
【理論】 経路積分法は、N 個の量子的な粒子を N×P 個の古典的な粒子(P：ビーズ数)として扱うことで、原子核の量子性を表現する。本研究では、配置生成に MD 法、配置の採択・棄却に MC法を用いるハイブリッドモンテカルロ法を採用した[6]。G-C pair に対して 150, 300, 450 K、A-T pairに対しては 150, 300 K の温度条件で計算した。ビーズ数は 150 K で P = 32、300 と 450 K では P =16、ステップ数は 400 000 である。また、核の量子効果の寄与を考察するために、核を一つの電荷として取り扱う古典シミュレーションも実行した。ビーズ数は P = 1、ステップ数は 1 920 000である。電子状態はすべての計算において PM6-DH+[8]で評価した。

【結果・考察】 図１の A-T pair と G-C pair の N1 原子と H1 原子の距離(RN1H1)を表１に示した。平
衡構造(Eq.)で の値とともに、各

Table 1 Average values and statistical errors of N1H1 distance in A-T and G-C pairs
together with the equilibrium values [Å]温度における古典シミュレーショ
ン(Cla.)と量子シミュレーション(Qm.)から計算された期待値を

Eq.
A-T pair
RN1H1
G-C pair
RN1H1

Cla.

150 K
Qm.

Cla.

300 K
Qm.

1.030

1.031(0) 1.040(0) 1.030(0) 1.040(0)1.029

1.030(0) 1.040(0) 1.030(0) 1.039(0)

Cla.

450 K
Qm.

1.031(0) 1.039(0)示した。Eq.での値と Cla.の期待値はほとんど変わらないのに対して、核の量子効果を考慮した Qm.
の期待値は大きい。これは、RN1H1 に対する非調和ポテンシャル曲線において、核の量子効果が反映されたためである。また、150 K の低温だけでなく 450 K においても核の量子効果は NH 距離に

Distribution (arbitrary unit)

影響を及ぼすことがわかる。
(a)

(b)

δH1

Cla.
Qm.

(c)

δH2

[Å]

δH3

Fig.2 One-dimensional distribution of (a)δH1, (b) δH2 and (c) δH3 on G-C pair at 300 K, respectively

300 K の G-C pair における水素結合に関わるプロトンの相対的な位置(δHi = RNiHi－RHiOi/HiNi’,
i=1-3)の一次元分布を図２に示した。実線が古典シミュレーション、破線が量子シミュレーションの分布を表している。図２において Cla.と Qm.を比べると、Qm.はδHi = －0.5 から 0 近傍に分布をもつ
ことがわかる。すなわち、核の量子効果を考慮するとプロトンは水素結合の中央に位置しやすくなる。
また、その傾向は H1 よりも H2 と H3 において強いこともわかる。この結果は、孤立気相中における塩基対間の二重プロトン移動が協奏的な H2 と H3 の反応であるというモデル[9]を支持する。
次に、150 と 300 K における古典シミュレーションと量子シミュレーションから得られたトラジェクトリ
に対して主成分分析を適用した。主成分分析では、各原子の変位から分子全体の変位を表す新しい合成変数を求めることができる。A-T pair のトラジ
ェクトリから得られる四つの主成分の割合を表２に示した。その四つの主成分は、A-T pair の振動数が最も 低 い 四 つ の 振 動 モ ー ド (buckle, propeller,

Table 2 The ratios of the contribution of thefirst four principal components in A-T pair at
150 and 300 K from classical and quantum
simulations. Unit in %.

opening, stagger)[10]に対応する。各シミュレーショ

description

ンにおいて分子運動の約６割を占める buckle とA-T pairbucklepropeller に注目すると、Cla.に比べて Qm.ではpropeller

buckle の割合が減り、propeller の割合が増えること

openingがわかる。これより、核の量子効果は propeller の分

stagger

150 K
Cla. Qm.

300 K
Cla. Qm.

39
21
10
4

38
24
11
4

33
25
1
3

20
34
10
5

子振動を促進させるということが考えられる。
プロトンの相対的な位置に関する温度依存性の解析や周囲環境を考慮した結果は当日ポスター
にて発表する。
[1]J. Sutor, Nature 195, 68-69 (1962) [2]J. Watson and H. C. Crick, Nature 171, 737-738 (1953) [3]S. Urashima, et al. J.
Phys. Chem. A. 114, 11231-11237 (2010) [4]J. Sponer, et al. J. Am. Chem. Soc. 126, 10142-10151 (2004) [5]K. Yagi, et al.
Chem. Phys. Chem. 10, 1442-1444 (2009) [6]K. Suzuki, et al. J. Chem. Phys. 132, 114108 (2010) [7]M. Daido, et al. Theor.
Chem. Acc. 130, 385-391 (2012) [8]M. Korth, J. Chem. Theory Comput. 6, 3808-3816 (2010) [9]J. P. Ceron-Carrasco, et al. J.
Phys. Chem. A. 113, 10549-10556 (2009) [10]A. R. Srinivasan, et al. Biophys. Rev. 1, 13-20 (2009)

4P105
ABINIT-MP(X) による京での FMO 計算（東大生産研 1 , 立教大理 2 , みずほ情報総研 3 , 産総研 4 , お茶大院人文創成科学 5 , 神戸大院シス
情 6 ）沖山 佳生1 , 望月 祐志 1,2 , 福澤 薫 1,3 , 古明地 勇人 4 , 森 寛敏 5 , 渡邉 千鶴 1 , 田中 成典 6

FMO calculations with ABINIT-MP(X) on the K computer(Univ. Tokyo, IIS1 , Rikkyo Univ.2 , Mizuho-IR3 , AIST4 , Ochanomizu Univ.5 , Kobe Univ.6 )Y. Okiyama1 , Y. Mochizuki1,2 , K. Fukuzawa1,3 , Y. Komeiji4 , H. Mori5 , C. Watanabe1 , S. Tanaka6【序論】フラグメント分子軌道法 (FMO 法)[1] 計算プログラムである ABINIT-MP(X)[2,3] を用いて、これまで様々なタンパク質をはじめとする大規模生体分子系の計算が行われている。

Fortran 言語および標準的な並列化手法を用いて開発されて来た ABINIT-MP(X) は汎用性のあるプログラムとして PC クラスタから個人のワークステーションに至る様々な計算環境で動作し、特に地球シミュレータ (ES2) においては、ベクトル機向けチューニングによって、世界最大規模の量子化学計算を非常に高効率で実行している (インフルエンザウイルス膜タンパク質とその抗体複合体 36160 原子系に対する FMO2-MP3/6-31G が 1024VPU を用いて 5.8 時間で完走[4])。2012 年秋から本格稼働されるペタスケールスパコンである京においても、超並列性とハードウェアの特徴を活かしてプログラムの性能を発揮し、実用性のある大規模 FMO 計算を実行することが重要である。
その実用上の一つの計算モデルが、高次多体展開 FMO 計算による高精度詳細相互作用解析である。我々のグループでは 4 体補正まで取り込んだ FMO4 計算 [5]FMO4
Etotal
=

∑
I

EI0 +

∑
I>J

∆ẼIJ +∑
I>J>K

∆ẼIJK +∑∆ẼIJKL(1)

I>J>K>Lを開発し、従来より細かいフラグメント分割手法と組み合わせることで、精度を落とすことなく、さらに詳細なフラグメント間相互作用 (IFIE) 解析の提案を行っている。このような詳細な分割が機能部位ごとの相互作用解析を可能とし、近年注目されている創薬手法の一つである

fragment-based drug design (FBDD) の観点からも、FMO 法の創薬への応用上重要な役割を果たすことが期待される。その一方で、多体展開項を 3 体、4 体へと拡張することで、SCF なら
びに高次相関計算の対象となるフラグメントペアの組み合わせが増大し、系全体での計算コストが問題になってくる。しかし、フラグメントサイズの縮小とそのようなタスクリストの増加は京スパコンのような省メモリ型超並列計算機には上手くフィットし、2 体計算まででは満たすことのできなかった多数の並列ノードを利用することで、膨大なフラグメントペアを一気に処理することが可能となる。
また、超並列性を活かした計算モデルとして、分子動力学 (MD) と組み合わせ、多数の構造サンプリングを実行することが考えられる (FMO-MD [6])。複数のジョブを並列的に処理することによって、これまで取得が困難であった量のアンサンブルを用いた統計的な議論が可能になれば、大規模分子系における一つの大きなブレイクスルーを迎えることが期待される。

本発表では、京の試験利用期間のトライアル枠 (192 ノード／ 1536 コア、1 時間) で行ったテスト計算の結果について報告を行う。
【結果と考察】ここでは京を用いた FMO4 計算の実行結果を紹介する。テスト計算には半径 15Åで発生させた droplet 構造の水クラスタ (分子数 418) を用いた (図 1)。また MP2 レベルの計算
を基底関数 6-31G と 6-31G**について行い、フラット MPI による並列化、フラグメントごとに割り当てたコア数を 1 とした。ここで、それぞれ基底関数の数が 5434 および 10450 であり、分極関数の追加による積分コスト増を考慮しても、問題の規模が約 2 倍の増加に対して、計算時
間も 2 倍強と、ほぼ線形にスケールしていることが分かる (表 1)。また、両基底で SCF の方がMP2 よりもコストが高いことも見てとれる。その他の応用事例やハイブリッド並列による実行
性能は当日ポスターセッションで報告する。
表 1: FMO4-MP2 計算の計測時間とタスクリスト数

6-31G

6-31G**

Task list(NF Cnmers )Monomer SCF

60.5

85.1

418

(418)

Monomer MP2

3.3

3.3

Dimer SCF

13.4

19.3

5219

(87153)

Dimer MP2

11.5

11.4

Trimer SCF

65.1

123.3

21180

(12085216)

Trimer MP2

30.2

64.8

Tetramer SCF

216.4

463.1

40811

(1253841160)

Tetramer MP2

117.3

340.9

FMO (Total)

600.9

1301.0【 謝 辞 】本 研 究 は 東 大 生 産 研 RISS プ ロ ジ ェ ク ト の 支 援 の も と で 行 わ れ て い る 。
【参考文献】[1] K. Kitaura, et al., Chem.

Phys.

Lett., 312 (1999) 319.

[2] T.

Nakano, et al., Chem. Phys. Lett. 351(2002) 475.

[3] Y. Mochizuki, et al.,Chem. Phys. Lett.

457 (2008) 396.

[4] Y. Mochizuki, et al., Chem. Phys.
Lett. 493 (2010) 346. [5] T. Nakano, etal., Chem. Phys. Lett. 523 (2012) 128.
[6] Y. Komeiji, et al., InTech Open Access Review Journal hhttp://bit.ly/
HkefpIi.
図 1: 水クラスタ (H2 O)418 の構造

4P106

Grb2 SH2 ドメイン-阻害剤間の結合親和性のFMO4 法と新規フラグメント分割法に基づく CH/π相互作用解析（東大・生産研 1，みずほ情報総研 2，国立衛生研 3，立教大・理 4，微生物化学研究所 5，CHPI
研究所 6）渡邉 千鶴 1，福澤 薫 1, 2，中野 達也 1, 3，望月 祐志 1, 4，梅沢 洋二 5，西尾 元宏 6

CH/π interaction analysis
for binding affinity between Grb2 SH2 domain and inhibitor
based on FMO4 method and new fragmentation scheme（Univ. Tokay IIS1，Mizuho-IR2，NIHS3，Rikkyo Univ.4，IMC5，CHPI Institute6）
Chiduru Watanabe1， Kaori Fukuzawa 1, 2, Tatsuya Nakano 1, 3， Yuji Mochizuki 1, 4，
Yoji Umezawa5，Motohiro Nishio 6【序】Growth receptor binding protein 2 (Grb2)の Src Homology 2 (SH2)ドメインは、成長因子受容体のチロシンリン酸化された配列に結合するタンパク質である。その阻害剤である
n 員環構造（図１(a)：pTyr+1 部分, n=3～7）を持つトリペプチドは、n 員環のサイズに依存して Grb2 SH2 ドメインとの結合親和性が変化することが確認されている[1]。その要因として、n 員環と Grb2 SH2 ドメインの Phe108 の側鎖の芳香環との CH/π相互作用（主に分散力が重要）などの相互作用が指摘されている。しかし、それらの指摘は X 線結晶構造からの知見であり、量子化学計算に基づいて分散力を見積もるような解析はされていない。本研究では、最近開発された 4 体補正フラグメント分子軌道（FMO4）法[2]と新規フラグメント分割法[3]を用いて、Grb2 SH ドメインと各阻害剤との相互作用を官能基単位のフラグメント間相
互作用エネルギーにより見積もり、それらの値を実験値と比較することで n 員環のサイズに依存した結合親和性の変化の原因を調べた。
【理論】本研究では、ABINIT-MPを用いてFMO4計算を行った[4]。Grb2 SH2ドメイン-阻害剤複合体を用いて、阻害剤を３分割（図１(ａ)）、タンパク質は主鎖-側鎖分割（図１(ｂ)）
した場合のFMO4計算(MP2/6-31G)を行い、4体項を取り入れたフラグメント間相互作用エネ
ルギー(FMO4-IFIE)を解析した。FMO4-IFIEを以下の式で定義する。

1
1
~
~
~
ΔE IJFMO4 = ΔE IJ + ∑ ΔE IJK + ∑ ΔE IJKL
6 K >L
3 K
ここで、I, J, K, Lはフラグ
メントの番号を表す。2体項（ａ） 阻害剤分割（ｂ） 主鎖‐側鎖分割までを計算するものを
FMO2-IFIE（通常のFMO
計算におけるIFIE）、3体項まで計算するものを
FMO3-IFIEと呼ぶ。また、

pTyr

pTyr+1 pTyr+2

FMO4-IFIE解析による
図１:新規フラグメント分割方法

CH/π相互作用の見積もりに対して検証を行

‐5.00
1.0
‐6.00

なうため、CHPI解析[5]、軌道エネルギー解析（FILM）[6]も併せて行った。
【結果と考察】各阻害剤（n=3～7）の pTyr+1
フラグメントと Grb2 SH2 ドメインとのFMO4-IFIE の和を求め、実験値(log Ka)との相関を調べたところ、IFIE の分散力の項

FMO4‐IFIE‐corr sum

‐7.00

10.0

100.0

三員環

‐8.00

四員環

‐9.00

五員環

‐10.00‐11.00‐12.00六員環‐13.00七員環‐14.00log Ka（FMO4-IFIEcorr）と実験値の間に良い相関があることが確認された（図２）。さらに、

R= ‐ 0.80

図２: FMO4-IFIE 和と実験値 Ka の相関

Grb2 SH2 ドメインと阻害剤（六員環）間のCHPI 解析の結果を図３に示す。特に、阻害剤の n 員環と Grb2 SH2 ドメインの Phe108との間にみられる CH/π水素結合の数は n 員

Lys109環のサイズに従って増加し、Arg67, Lys109,
Trp121 との間の CH/π水素結合の数は n 員環のサイズに依存しないことが確認された。
CHPI 解析の結果は、pTyr+1 フラグメントと各アミノ酸残基間の FMO4-IFIE の結果とも

Trp121

一致し、結合親和性の変化の主な原因が n 員
環と Phe108 との CH/π相互作用であること

Arg67

が特定された。発表当日には、FILM 解析による検証も報告する。

Phe108

図３: Grb2 SH ドメインと阻害剤間の CHPI 解析【謝辞】
本研究は、文部科学省次世代 IT 基盤構築のための研究開発「イノベーション基盤シミュレー
ションソフトウェアの研究開発」プロジェクトにおいて実施された。
【参考文献】[1] J. M. Myslinsk, J. E. DeLorbe, J. H. Clements, and S. F. Martin, J. Am. Chem. Soc.,
133, 18518–18521 (2011)[2] T. Nakano, Y. Mochizuki, K. Yamashita, C. Watanabe, K. Fukuzawa, K. Segawa, Y.
Okiyama, T. Tsukamoto, S. Tanaka, Chem. Phys. Lett., 523, 128–133 (2012).
[3] C. Watanabe, K. Fukuzawa, Y. Okiyama, Y. Mochizuki, T. Nakano et. al., to besubmitted.
[4] BioStation 6.0: ABINIT-MP and BioStation Viewer. The program package is due tobe released.
[5] Y. Umezawa, M. Nishio, Bioorg. Med. Chem., 6, 493-504 (1998).
[6] T. Ishikawa, Y. Mochizuki, S. Amari, T. Nakano, H. Tokiwa, S. Tanaka, K. Tanaka,

Theor. Chem. Acc., 118, 937-945 (2007).

4P-107
一酸化窒素還元酵素 cNOR 及び qNOR における一酸化窒素還元反応機構の理論的研究（京大福井謙一セ 1、エモリー大化学 2）
畑中 美穂 1、阿部 恵理花 1、諸熊 奎治 1,2

Theoretical study of the mechanism of NO reduction
in bacterial nitric oxide reductase cNOR and qNOR
(Fukui Institute for Fundamental Chemistry Kyoto University1, Emory University2)
Miho Hatanaka1, Erika Abe1, Keiji Morokuma1,2【序】
嫌気性呼吸のプロセスの中で、有毒な NO を無毒な N 2 O に変換する酵素を一酸化窒素還元酵素(NOR)と言う。この酵素は、好気性呼吸において酸素を水に還元する酵素(CcO)と同じ祖先から進化したため、嫌気性呼吸から好気性呼吸への進化の謎を解き明かす鍵を握ると考えられている。NOR から CcO への進化の過程で大きく変化した点は二点ある。一つ目は、反応中心に含まれる金属が、NOR ではヘム鉄・非ヘム鉄であったのに対し、CcO ではヘム鉄・非ヘム銅に変化したことである。二つ目は、NOR では電子とプロトンの供給源が、共に細胞の外側にあったのに対し、CcO では、プロトンの供給源だけが、細胞の内部に変わったことにより、プロトンポンプの機能が生まれたという点である。
しかし、最近になって、これまでに結晶構造が明らかにされていた NOR の一種である cNORとは別の qNOR の結晶構造が報告され、qNOR の全体構造は、cNOR や CcO と似ており、反応中心も cNOR と似ているものの、プロトンの通り道が cNOR とは異なり、細胞の内側からのびていること、その通り道は CcO と似ていることが明らかにされた。このことから、qNORには、CcO におけるプロトンポンプの原型ともいえる“試作品”が存在しており、cNOR
よりも、より CcO に近い酵素であることが推測された。
cNOR と qNOR の違いは、プロトンの通り
道だけではない。両者の反応活性点は、同じ
金属、アミノ酸残基で形成されているものの、一つのアミノ酸残基(Glu)が cNOR では非ヘム鉄に配位し、図 1 の Cmp1 の段階で非ヘム
鉄が五配位になっているのに対し、qNOR では、Glu が非ヘム鉄から大きく離れた場所に位置するため、非ヘム鉄が(水の配位を除く

図 1：提唱されてきた反応機構と)三配位の状態になっている。そのため、
qNOR では、NO が非ヘム鉄に配位する水分子と配位交換を起こせば、cNOR よりも多くのNO が非ヘム鉄に配位し得るため、異なる反応機構を経る可能性がある。
従来、ヘム鉄・非ヘム鉄を含む反応中心で NO 二分子が N2O へ還元される機構として、cis:b3
機構、cis:FeB 機構、trans 機構などが提唱されてきた。cNOR のモデル錯体を用いた DFT 法による先行研究 1,2 では、反応の初期段階で非ヘム鉄と NO 二分子による五員環構造(図 1 の Cmp7)

を形成するという段階までは一致するものの、その後に続く N-O 結合の解離には、異なる機構が提唱されていた。qNOR の場合は、非ヘム鉄の配位数の違いから、cNOR では不利な機構と考えられている cis:FeB 機構も十分起こり得ると予測できるが、これまで検討されていなかった。そこで我々は、モデル錯体を用いた DFT 法、タンパク全体を考慮した ONIOM 法及び、反応経路を自動的に探索するために人工力誘起反応(AFIR)法を用いることで、cNOR、qNORの NO 還元機構違いを調べ、非ヘム鉄の配位数や、反応活性点の周囲のタンパクが、反応機構にどのように影響を及ぼすのかを明らかにすることを目指す。
【計算方法】
計算モデルを結晶構造(cNOR は PDB ID =3O0R、qNOR は PDB ID =3AYF)を元に作成し、活性中心を取り出したモデル錯体に対するDFT 計算と、酵素全体に対する ONIOM (DFT : MM)計算を行った。
DFT 計算には、UB3LYP/[LANL2DZ(Fe),6-31G*]を用い、図２の*の原子を固定して構造最適化を行った。また、ONIOM 法の MM 計算には Amber 力場を使用した。
図 2：モデル錯体【結果】

まず、cNOR における NO 還元機構について DFT 法を用いて調べたところ、図 3 のように
非ヘム鉄と NO 二分子が五員環構造を形成した後(cis:b3 機構)、一プロトン、一電子の N2O2への移動が起き、N-O 結合が開裂することで、N2O が生成するという機構が、最も反応障壁の実験値をよく再現することを明らかにした。この反応の中で律速となるのが、N-O 結合の開裂である。非ヘム鉄と NO 二分子の五員環（図１Cmp7）から直接 N-O 結合が開裂する場合、その反応障壁は 31 kcal/mol もあり、実験値 13.3 kcal/mol と比べると明らかに高い値であった
が、電子、プロトンを一つずつ付加した場合の N-O 結合開裂の反応障壁が 11.0 kcal/mol と、
実験値と良い一致を示すことが分かった。また、ONIOM 法を用いて周囲のタンパクを考慮しても、N-O 開裂の反応障壁に大きな変化はないことも分かった。

図 3：cNOR における NO 還元機構
さて、
cNOR では、
非ヘム鉄に配位する Glu がプロトンを N2O2 に渡す役目を果たすが、qNOR
では、この Glu が非ヘム鉄から遠くに位置しているため、プロトンの受け渡しができないと考えられる。そこで、改めて AFIR 法を用い、qNOR のモデル錯体に NO 二分子を付加させる
人工力を加えて、反応経路を探索したところ、cNOR では得られなかった反応経路を複数得た。発表では、cNOR と qNOR の反応機構の違いについて議論する。
【参考文献】[1] M.R.A. Blomberg, et al, Biochemistry, 2012, 51, 5173. [2] L.M. Blomberg, et al,
Biochim. Biophys. Acta 2006, 1747, 240.

4P108

分子動力学計算における相互作用計算の高速化手法の評価（金大院・自然科学） 吉田津日門、齋藤大明、川口一朋、長尾秀実

Evaluation of acceleration method of interaction
calculation in molecular dynamics calculation(Kanazawa Univ.) Tsuhito Yoshida, Hiroaki Saito, Kazutomo Kawaguchi, Hidemi Nagao

1 [序]分子動力学 (MD) 計算は、タンパク質などの生体分子系でのシミュレーションに使われている。しかし、MD 計算は計算量が粒子数の二乗に比例するため、生体分子系の様な大規模系では計算時間が膨大になる。それを解決するためには、計算機の性能を上げる他に、計算プログラムを見直すことでも高速化が可能である。そこで、MD 計算の計算時間の大部分を占めている LJ 相互作用とクーロン相互作用の計算の高速化を行い、それぞれ複数ある高速化手法について計算時間を比較する。本研究では水分子モデルの一つである
TIP3P を用いる。TIP3P を用いることで LJ 相互作用、クーロン相互作用の計算ができ、また粒子数を増やすことでタンパク質系の計算と同様な計算量で計算の実行ができる。その結果からタンパク質系のシミュレーションにおいてどのような手法が有効であるか考察する。
2 [方法]
TIP3P を用いて粒子数、カットオフ距離を変化させた場合において各高速化手法の相互作用計算時間を比較する。計算条件は周期境界条件下で密度一定の基本セル中に水分子があるとして、温度は 300K とする。カットオフ距離を変える場合は分子数 256、粒子数を変える場合はカットオフ距離 10Å で計算する。LJ 相互作用の様な短い距離で収束する相互作用計算では、最初にカットオフを行う。しかし、全粒子対についてカットオフ処理をすると、相互作用が明らかに働かない粒子対の距離も求めなければならない。以下にその距離計算の回数を減らす方法を示す。LJ 相互作用の高速化手法は 2 つあり、一定距離内の粒子を数ステップごとに登録する帳簿法と、セルを区切って計算範囲を減らすセルイ
ンデックス法がある。
2.1 帳簿法
まず、予めカットオフ距離 (Rc) よりも広い距離 (Rs) 内にいる粒子を登録する。そして登録した粒子に対してのみカットオフ距離と比較する。粒子の登録は数ステップごとに更新をする。粒子登録の更新時間の間隔（∆t ）と Rs の大まかな設定は次の式 (1) で求める
[1]。
∆t ≤

Rs − Rc
2VM AX(1)

VM AX は粒子の最大速度を表す。
2.2 セルインデックス法
まず、セルの各軸方向にカットオフ（あるいはそれ以下）の距離で小さいセル（サブセル）に区切る。次に各サブセル内にある粒子の数と座標を登録する。そして対象の粒子が登録されているサブセルと隣り合ったサブセル内の粒子との距離を計算し、カットオフ距

離との比較をする。この方法はカットオフ距離がセルの長さの 4 分の 1 以上の場合、すべてのサブセルについて粒子間距離を求めなければならなくなるため、カットオフ距離がセルの長さの 4 分の 1 以下でなければ効果がない。
2.3 帳簿法・セルインデックス法両方を用いた方法
まず、帳簿法で用いるカットオフ距離よりも広い距離でサブセルに区切り、セルイン
デックス法と同様の操作を行う。次に、それぞれ区切ったサブセルにおいて、帳簿法を用いて粒子の登録をする。その後は毎ステップ、サブセル内の粒子を更新し、式 (1) で求めた時間 ∆t ごとに帳簿の更新をする。
2.4 PME(Particle Mesh Ewald) 法
PME 法はクーロン相互作用の高速化手法の一つである。周期境界条件でのクーロン相互作用を計算するには、エワルド法を用いる必要がある。エワルド法とは、セル内の粒子対は実空間の計算をし、同セルでない粒子対は逆格子空間の計算を行う方法である。そこ
でさらに逆格子空間の計算を高速フーリエ変換 (FFT) を用いて計算する方法が PME 法である。

3 [結果と考察]
LJ ポテンシャルについての計算結果を示す。カットオフ距離と粒子数について高速化手法を導入していない場合と、帳簿法（以下、手法１）、セルインデックス法（手法２）、両方を用いた場合（手法３）で計算した図をそれぞれ図 1,2 に示した。カットオフ距離を変化させた場合と粒子数を変化させた場合でどの方法も変更前と比べて高速化できている
ことが分かる。図 1 から、カットオフ距離が大きくなると手法１の方が他の手法より速くなり、さらに大きくなるに連れて手法 1 と手法 2,3 との差が大きくなっている。したがってカットオフ距離がセルに対して大きく取らなければならない場合、手法１を用いると良いことが分かる。また、図 2 から、粒子数が多い場合は手法 2 を用いると良いことが分かる。クーロンポテンシャルについての高速化手法との比較は当日会場で発表する。
2.0e−03

0.06
Book keeping
Cell index
Book+CellNone

1.8e−03

Book keeping
Cell index
Book+Cell
None
0.05

1.6e−031.4e−03
0.04

time[sec]time[sec]

1.2e−031.0e−03

0.03

8.0e−04
0.02
6.0e−044.0e−040.01
2.0e−04
0.0e+00

0
5

10

15

20
cutoff range[Å]

25

30

35

0

500

1000

1500
Particle number

2000

2500

3000

図 1: カットオフ距離に対する１ステップあ 図 2: 粒子数に対する１ステップあたりの計
たりの計算時間
算時間

4 [参考文献]［1］ H.Watanabe, M.suzuki, and N.Ito Eﬃcient Implementions of Molecular Dynamics
Simulations for Lennard-Joned Systems Progress of Theoretial Phisics Vol.126 August
2011 203-235

4P-109
α ヘリックスにおける L-と D-アミノ酸残基の双極子モーメントに関する量子化学的研究（広島大院理、広島大 QuLiS）三枝俊亮, 相田美砂子

Quantum chemical study on the dipole momentsof L- and D-amino acid residues in α-helix(QuLiS and Graduate School of Science, Hiroshima University) Shunsuke Mieda,
Misako Aida
【序】グリシンを除くアミノ酸のα炭素はステレオジェニック中心であるため、L 体
および D 体をとりうる（図１に Serine の L 体および D 体を示した）。自然界に存在するほぼ全てのアミノ酸は L 体をとるが、どのような化学的性質のために L 体が優位になるのかは、これまで明らかにされていない。我々は以前の研究で、タンパク質中に存在するアミノ酸残基の双極子モーメントを計算し、多くの L－アミノ酸残基の双極子モーメントが D－アミノ酸残基の双極子モーメントに比べ大きいことを見いだした。すなわち、双極子モーメントの大きさの違いが、L 体の優位性の一因となった可能性がある。L 体と D 体では側鎖側鎖の出る方向が異なっており（図１）、双極子モーメントの大きさの違いはこのためであると考えた。双極子モーメントの違いが側鎖部の違いによるということを示すた
めに、アミノ酸残基の双極子モーメントを主鎖部と側鎖部に分割する方法を提案する。
【計算・解析方法】アミノ酸残基の双極子モーメントを原子電荷と座標の積の総和で表わす。双極子モーメントの計算に用いる原子電荷は、MP2/6-311G**レベルで計算
し、Natural Population Analysis を用いて求めた。我々は図２に示したような、α炭素の原子電荷を２つに分割する手法を提案する。

図２

アミノ酸残基の原子電荷を主鎖部と側鎖部に分ける。その時に主鎖部と側鎖部を単純に分割するのではなく、原子電荷の総和が各部でそれぞれ 0 になるようにα炭素の原子電荷 qcを qc1 および qc2 に分割する。なお、分割したα炭素の原子電荷は、いずれもα炭素の位置（rcに置く。このように、主鎖部の原子電荷の総和と、側鎖部の原子電荷の総和がいずれも０とすることで、双極子モーメントの計算時に原点の位置を考慮する必要がなくなる。分割した原子電荷とそれぞれの原子の座標を用いて双極モーメントの計算を行った。
図３に alanine を例として、具体的な分割方法を示す。まず、alanine の原子電荷を計算する（図３(a)）。次に、α炭素の原子電荷を分割し（図３(b)）、分割したα炭
素の原子電荷を用いて各部の双極子モーメントを計算する（図３(c)）
【結果と考察】図３(c)に L-alanine と D-alanine の双極子モーメントを示す。L 体の双極子モーメントのほうが D 体の双極子モーメントよりも約 0.28Debye 大きい。双極子モーメントを主鎖部と側鎖部に分割したところ、主鎖部の双極子モーメントの大きさは L 体も D 体もほぼ等しい事がわかった。一方、側鎖の Y 軸方向の双極子モーメントは L 体が 0.198 Debye、D 体が－0.145 Debye であり、これらの双極子モーメ
ントの差が全体としての双極子モーメントの違いを生じさせる。また、valine についても同様に双極子モーメントの分割を行った結果、主鎖部の双極子モーメントはalanine の場合とほぼ同じ値であることがわかった。すなわち、主鎖の双極子モーメ
ントは同じ二次構造をとっている限り、側鎖の種類にかかわらずほぼ同じであり、光学異性体およびアミノ酸残基の種類による違いは、側鎖によるものである。

図３

精度記録型量子化学計算

4P110（北教大・釧路）小原繁Accuracy-Recording Quantum-Chemical Calculations(Hokkaido University of Education, Kushiro Campus)

1

Shigeru Obara【序】

まうことが分った。別の言い方をすると、より正確な精コンピュータの発達に伴い量子化学計算の研究対象 度算出には「精度が増す要因」を探す必要があることは大型化して来ている。その一方で計算そのものは依 が分った。
然として倍精度実数を使って行われており、精度不足「確実に精度が増す要因」を四則演算だけから探しの発生が懸念される。精度不足発生の有無そのものは、 出すことはできない。考え得るものは「確率的に精度
別途精度の高い計算を行って比較することにより確認 が増す要因」だけである。この要因として、本研究ではすることはできる。しかし、このような別途計算を常 以下の三つの演算を取り上げ「精度追加」を行った。
に行うことは不便である。それに、きめ細かな解析を
1. 減算精度落ち時追加 減算精度落ちにより n ビッ
行うことも容易でない。通常の量子化学計算を行いなト残ったとする。この場合、精度は n ビットへ減っ
がら計算精度も常に記録していくことができるとこれたと考えるのではなく、さらに下位の 1 ビット（0らの不便が解消する。この観点から、どのように精度か 1 かのいずれか）が確率的にある程度正しいのを記録していけば良いのか等を検討して本研究室で開で、(n + α1 ) ビットの精度へ減ったと考える。この発した高精度量子化学計算プログラムに精度記録機能
α1 が「精度追加」である。
を追加した。これについて報告する。
2. 加算繰り上り時追加 例として有効ビット数が 32 【有効ビット数の増減要因】
ビットの下記 2 進数の加算を考える。
理系の数値計算では有効数字という概念が用いられ、
101
その桁数不足発生の有無が常に注視される。この場合の有効数字は 10 進数のことであるが、コンピュータで
(1)
+ 100
は 2 進数が使用されるので「2 進数の有効数字」という1001
意味で「有効ビット」と表記し、また、「その桁数」という意味で「有効ビット数」と今後表記する。
和の値は繰り上りにより 4 ビットになっている。こ有効ビット数が減少する主要因は近接数同士の減算の数の有効ビット数は元々の 3 ビットよりも多くが行われることである。より正確には、絶対値が近接なったと考えてよいだろう。つまり、加算結果が繰する同符号二数の減算、あるいは、異符号二数の加算
り上る時は有効ビット数が α2 ビットだけ増すと考が行われることである。
える。

−

3. 丸め繰り上り時追加 例として計算結果が下記の値だったとして、これを小数第 1 位で丸める場合を考える。
111.1 ⇒ 1000
(2)

これらの演算が行われると、有効ビットの上位が削ら丸めた値は 4 ビットになる。この数の有効ビット数れて下位だけが残り有効ビット数が減少する（上図）。
が元々の 3 ビットよりも多くなると考える。つまこの減少要因以外に有効ビット数が減少する要因を、り、丸め時に繰り上る時は有効ビット数が α3 ビッ
通常、省られることはない。一方、有効ビット数が増
トだけ増すと考える。なお、この繰り上りの起き加する要因については、通常、全く考慮されることがる頻度はかなり低い。しかし、一応、有効ビットない。
数が増す要因の一つとして掲げておく。
計算精度を決める以上の要因、つまり、「減算だけによる精度減少」により (18 × 18) 対称行列の対角化計算 3 【計算方法】（16 倍精度計算）を行い、計算結果の各要素の有効ビッ
上記の「精度追加」α1 , α2 , α3 によって有効ビット数
ト数を算出した。また、この算出結果の精度を確認する
を適切に算出することができるかどうかを確認するたために、高精度計算（32 倍精度計算）を別途行って得めに、Rn 原子の量子化学計算を行ない、重なり積分のた結果から精度を算出した。両者を比較すると、極めHouseholder 法による三重対角化における有効ビットて残念なことに、差が極端に大きかった。実際には、20数の算出を行なった。
ビット強しか精度落ちしていないのに、40 ビット強の量子化学計算には本研究室において Java 言語を使用
精度落ちがあると算出され（表 1 の 2 行目と最終行）、して開発した量子化学計算プログラムを使用し、16 倍「減算だけによる精度減少」では精度を低く算出してし精度（仮数部 496 ビット）の計算を解析対象データと

表 1: Rn 原子 16 倍精度計算の重なり積分を House-

holder 法で三重対角化したときの最低精度要素の有効ビット数の算出値と実際の値
表 2:

α1 と α2 の種々の値の組み合わせによる Rn

α値†

主対角要素

副対角要素

原 子 16 倍 精 度 計 算 の 重 な り 積 分 三 重 対 角 化 後 の減 0.00, 加 0.00, 丸 0.00

450
459
467
472
476

451
461
468
473
477

最低精度要素の有効ビット数の算出値と実際の値。

450
456
461
464
466

451
457
462
464
466

いは、減算桁落ち時追加精度を大きくして加算繰り上

450
450
450
450
450

451
451
451
451
451

α値†

主対角要素

副対角要素

減 0.00, 加 0.00

472.54‡473.58‡

450
456
461
464
466

451
457
462
464
466

459
466
469
472
473

461
467
470
472
473

467
471
475
476
478

468
472
475
477
478

472
476
478
479
480

473
476
478
480
481

476
478
480
481
482

477
479
481
482
483

472.54‡473.58‡

減 0.25, 加 0.00, 丸 0.00
減 0.50, 加 0.00, 丸 0.00
減 0.75, 加 0.00, 丸 0.00
減 1.00, 加 0.00, 丸 0.00
減 0.00, 加 0.00, 丸 0.00
減 0.00, 加 0.25, 丸 0.00
減 0.00, 加 0.50, 丸 0.00
減 0.00, 加 0.75, 丸 0.00
減 0.00, 加 1.00, 丸 0.00
減 0.00, 加 0.00, 丸 0.00
減 0.00, 加 0.00, 丸 0.25
減 0.00, 加 0.00, 丸 0.50
減 0.00, 加 0.00, 丸 0.75
減 0.00, 加 0.00, 丸 1.00†

この列の 3 つの数は α1 , α2 , α3 の値を表わす‡

32 倍精度計算から算出した実際の有効ビット数

実際の値に近い [α1 , α2 ] の値は [0.25,0.75](472,472)、[0.25,1.00](473,473)、および、[0.75,0.00](472,473) だった。これは、減算桁落ち時追加精度を小さくして加算
繰り上がり時追加精度を大きくする（前者２つ）、ある
がり時追加精度をなしにする、と実際に近い値を算出できると解釈できるが、より確実なものにするには一層の比較・検討が必要だろう。

減 0.00, 加 0.25
減 0.00, 加 0.50
減 0.00, 加 0.75
減 0.00, 加 1.00
減 0.25, 加 0.00
減 0.25, 加 0.25
減 0.25, 加 0.50
減 0.25, 加 0.75

した。また、精確な値を得るために 32 倍精度（仮数部992 ビット）の計算を行ない、これを真値として用いた。
なお、このとき数値計算上の誤差だけを解析に使用するため、解析対象データである 16 倍精度の重なり積分の各要素の仮数部下位に 496 ビットの 0 を追加することにより得た行列を 32 倍精度計算で用いる 32 倍精度
重なり積分とした。

4【精度追加と有効ビット数】

それぞれの「精度追加」α1 , α2 , α3 の種々の値と、それを使用して算出された重なり積分三重対角化行列の最低精度の主対角要素と副対角要素の有効ビット数を
表 1 に掲げる。
この表の最左列は α1 , α2 , α3 の値を表わす。値が 0.00
は追加する有効ビット数が零ビットであることを表わ
し、1.00 は追加する有効ビット数が 1 ビットであることを表わす。表の中央列と最右列はそれぞれ主対角要素と副対角要素の中で最低精度の有効ビット数の算出値である。また、最終行には 32 倍精度計算から得られた当該要素の実際の有効ビット数を掲げる。
表から α1 や α2 の値によって算出された有効ビット数が大きく変わることがわかる。一方、α3 では当初の予想通り算出有効ビット数に影響が現れない。
α1 と α2 の種々の値の組み合わせによって三重対角化
行列の最低精度要素の有効ビット数の算出値がどう変わるをまとめたのが表 2 である。詳細は当日発表する。

減 0.25, 加 1.00
減 0.50, 加 0.00
減 0.50, 加 0.25
減 0.50, 加 0.50
減 0.50, 加 0.75
減 0.50, 加 1.00
減 0.75, 加 0.00
減 0.75, 加 0.25
減 0.75, 加 0.50
減 0.75, 加 0.75
減 0.75, 加 1.00
減 1.00, 加 0.00
減 1.00, 加 0.25
減 1.00, 加 0.50
減 1.00, 加 0.75
減 1.00, 加 1.00†

この列の 2 つの数は α1 と α2 の値を表わす‡

32 倍精度計算から算出した実際の有効ビット数

4P-111
SnPc, PbPc 及び Sn(Pc)2 の電子構造と励起エネルギーに関する理論的研究（山口大院理工 1，熊本大院自然 2，横市大院生命ナノ 3）
隅本倫徳1，本多輝行 2，川島雪生 3，堀憲次 1，藤本斉 2

Theoretical investigation on the electronic structures and excitation energies of
MPc and M(Pc)2 (M = Sn and Pb)(Yamaguchi Univ.1, Kumamoto Univ.2, Yokohama City Univ.3)
Michinori Sumimoto1, Teruyuki Honda2, Yukio Kawashima3, Kenji Hori1, Hitoshi Fujimoto2【序】N

フタロシアニン(Pc)類は、クロロフィルやヘモグロビンと構造が似ていることから生態系のモデルとして注目され、また、化学的、熱的にも安定であることから機能性色素として古くから利用されてきた。これまで、金属フタロシアニンとしていろいろな化合物が合成されており、中

N

N
M

N

N
N

N
N

心金属を変えることにより、その性質が変化することが知られている。
通常、金属フタロシアニンは平面構造をしており、D4h の分子対称性を持

Scheme 1.

っている (Scheme 1)。しかしながら、Pc 環の中心に入りきらないよ
うなイオン半径の大きい Pb や Sn が配位した場合、バトミントンの羽根のようなシャトルコック型構造になることが知られている(Scheme 2)。SnPc 及び PbPc では、分子構造が似ているにもかかわらず吸収スペクトルに違いが見られる。このような物性の違いについて、詳細な理由は解明されていない。
本研究では、MPc と M(Pc)2 (M = Sn, Pb) の分子構造、電子構造および励起状態、また SnPc の結晶構造をモデル化した二量体から結晶構造の物性について詳細な理論的解釈を与えることを目的とした。
Scheme 2.
【計算方法】
計算は Gaussian 03 及び 09 プログラムを使用した。構造最適化、及びエネルギー計算には DFT法、励起エネルギー計算には Time-dependent DFT (TDDFT) 法を用いて計算した。TDDFT 法では、
M06、B3LYP、BPW91、BHandHLYP の汎関数を用いた。二量体の構造最適化及び励起エネルギー計算には、過去の研究 1,2 で良好な結果が得られた M06 汎関数を使用した。構造最適化及び励起エ
ネルギー計算には、Sn 及び Pb に LanL2DZ(d,p)、C, N, H に 6-311G(d)の基底関数を用いた。

【結果と考察】SnPc の構造最適化を行ったところ、C4v 対称構造が得られた。イオン半径の大きな原子がフタロシアニン環の中心に配位した場合、環内に入ることができずシャトルコック型構造を取る。また、イオン半径が大きくなれば、
フタロシアニン環からのとび出し方も大きくなる。これらの構造は以前に報告された X 線
結晶構造解析の結果とほぼ一致していた。得られた最適構造を用いて、励起エネルギー計算を行った。TDDFT 法で得られた励起エネル
ギーと実測した吸収スペクトルを Figure 1 に示す。SnPc は、1.7 eV 付近の吸収帯（Q バンド）が幅広くなっている以外は通常の D4h 対称構造を持った金属フタロシアニンと同様な吸Figure 1. Excited states of SnPc obtained by the

収スペクトルを示す。TDDFT 法を用いて得ら

TDDFT method.

れた励起エネルギーは、実測したスペクトルと良く一致している。3
次にＸ線結晶構造を参考にし、結晶構造中の二分子を取り出したよ
うな二量体構造の構造最適化を行った。SnPcの結晶構造には、飛び出した金属同士が内側を向くものと外側を向く
ものの二種類が考えら

Figure 2. Optimized geometries of two type SnPc dimers.

れる。得られた最適化構
造を Figure 2 に示した。それぞれ C2h 対称構造を持つこれらの二量体構造は、X 線結晶構造解析の結果と非常に似通った構造であった。得られた二量体構造を用いて、励起エネルギー計算を行った。これらの詳細な解析については当日の発表で報告する。
【文献】(1) M. Sumimoto, Y. Kawashima, D. Yokogawa, K. Hori, H. Fujimoto, J. Comput. Chem. 2011, 32, 3062.
(2) M. Sumimoto, Y. Kawashima, D. Yokogawa, K. Hori, H. Fujimoto, Int. J. Quantum Chem.,
DOI:10.1002/qua.24072.
(3) M. Sumimoto, T. Honda, Y. Kawashima, K. Hori, H. Fujimoto, Dalton Trans. 2012, 41, 7141.

4P-112Mean-field second-order Moller-Plesset
perturbation theory
Tran Nguyen Lan1 and Takeshi Yanai1, 2
1
The Graduate University for Advanced Studies
2
Institute for Molecular Science
We use the canonical transformation to include the dynamic correlation effect at theMP2 level on the molecular Hamiltonian under the one-particle (mean-field)
formalism. The cumulant decomposition plays the key role to do this purpose. By thisway, we will obtain a new correlated Hamiltonian. Therefore, the wave function andeigenvalues after relaxation due to correlation field can be obtained by resolving the
Schrodinger equation with this new Hamiltonian. Several tests have been done to show
advantages of our theory.

4P113
リチウムイオン電池正極における電解液分解反応に関する理論的研究（阪府大院理 1,RIMED2）櫻井耕司 1，麻田俊雄 1,2，小関史朗 1,2Theoretical study in decomposition of electrolyte on Li ion battery(Osaka Prefecture Univ1,RIMED2) Koji Sakurai1, Toshio Asada1,2, Shiro Koseki1,2【序論】
リチウムイオン電池は他の電池に比べ高容量、高電圧で、携帯電話やノートパソコンなどの電子機器に広く用いられており、負極に金属リチウムを含むため電解液としては非プロト
ン性の有機溶媒が用いられている。リチウムイオン電池の電解液のひとつである propylene
carbonate (PC) は電池の充電中に正極付近で二酸化炭素を発生して酸化的に分解されることが報告されており、電池の性能劣化の一因となっている。
そこで、本研究では量子化学計算と力場計算を組み合わせた quantum mechanical/molecularmechanical （QM/MM）molecular dynamics (MD) シミュレーションを行い、PC の分解反応の機構を理論的に解析し、分解されにくい新規分子の提案を行う。
(a)

(b)

E(kcal/mol)

3.0
TS1･+
16.7

1.4

TS4･+
PC1･+

P1･+

12.0

0

TS5･+
･+

PC･+
0

≈

e-

P2
-8.9

･+

P3
-35.4

e-

アセトン

プロパナール

図１ 真空中における PC の分解反応の反応機構(M06-2X/6-311++G(d,p))(a) 二酸化炭素分子の脱離過程, (b) 異性化反応【計算方法】
PC カチオンラジカル 1 分子の周りに溶媒として 670 個の中性 PC 分子と電解質 LiPF670 個を配置したものを基本セルとした。平衡構造を 5 ns 間の定圧 MM MD シミュレーションを用いて作成し、反応経路の最適化には QM/MM シミュレーションを行った。QM 領域には
M06-2X/6-311G(d) 法を、MM 領域には Amber 99 力場を用いた。また、自由エネルギー面上での反応経路の最適化には Free Energy Gradient (FEG) 法 1 と Nudged Elastic Band (NEB) 法
を組み合わせた FEG-NEB 法を用いた。FEG 法では QM 領域の自由エネルギー勾配∂A(rQM )/ ∂rQM は QM 原子に関するエネルギー勾配 ∂E(rQM , rMM )/ ∂rQM の MM 領域のアン
サンブル平均から次式に従って計算することができる。

∂A(rQM )∂E(rQM , 𝑟𝑀𝑀 )
=⟨
⟩
∂rQM
∂rQM

𝑟𝑀𝑀

なお、MD シミュレーションには Amber-gaussian interface1 を用いた。

【結果と考察】PC 分子の分解反応に関するモデル系として、すでに純粋な PC 溶媒に対する QM/MM シミュレーションにより得られた自由エネルギー面上の分解反応経路を解析してきた（図３）。そ・+の結果、溶媒効果により TS1 の反応障壁が 11.6 kcal/mol と真空中の値 16.7 kcal/mol と比べてかなり低くなる結果となった。これは反応の進行に伴い TS 付近で分子の双極子モーメントが大きくなり、溶媒からの安定化を強く受けることに起因していることを明らかとした。
一方、分解反応に溶媒中の電解質に由来する PF6-が直接影響する可能性について近年報告がなされた 3 。そこで PF6-イオンが分解反応に与える影響を解析した。PC・+に対して PF6-イ
オンが付加すると PC からプロトンが脱離した構造 2（図 4）が得られる。
（PC ＋ PF6- → PC(-H+) ＋ HF + PF5 ）PC(-H+)
PF5

PC・+から脱離したプロトン
図 3 PC 溶媒中におけるの自由エネルギーおよび双極子モーメントの変化

図 4 PCM//M06-2X/6-311++G(d,p)計算から得られた
脱プロトン化反応の生成物の三分子会合体

真空中および PC 溶媒中における脱プロトン化反応に伴う自由エネルギー変化を表１に示す。脱プロトン化反応により生成する三分子会合体は PC 溶媒中で 8.6 kcal/mol 自由エネルギ
ーが安定化する。また生成した三分子が解離した構造の自由エネルギーは－21.8 kcal/mol となり、後者のほうが自由エネルギー的に有利であることが明らかとなった。これは分子間の相互作用が弱くなることによりエントロピーが増加し自由エネルギーが低くなるためだと考えられる。
以上の結果をもとに三分子（PC(-H+), HF, PF5 ）が解離したモデルに基づいて電解液分解反応を解析し、純粋な PC 溶媒中の結果と比較検討した。詳細は当日発表する。
表１ 脱プロトン化反応に対する相対自由エネルギー
反応物（PC+ + PF6-）生成物(PC(-H+) + HF + PF5)
三分子会合体解離構造
真空中0.0-106.6-115.5
PC 溶媒中
0.0
-8.6
-21.8
（kcal/mol）【文献】
1. M. Nagaoka, N. Okuyama-Yoshida, and T. Yamabe, J.Phys.Chem. A,102,8202(1998).
2. 櫻井耕司ら, 第 15 回 理論化学討論会仙台 ポスター 2P25.
3. Xing, L.; Borodin, O.; Smith, G. D.; Li, W. J.Phys.Chem. A 115, 13896(2011).

4P114

有機ラジカル間の磁気的相互作用を利用した分子ワイヤの減衰定数 β の理論的計算(京大院工 1, 京大福井セ 2)

西澤 尚平 1, 長谷川淳也 1,2,

松田 建児 1

Theoretical investigation of the β value of the phenylene andphenylene ethynylene units by evaluating exchangeinteraction between organic radicals(Kyoto Univ.) Shohei Nishizawa, Jun-ya Hasegawa, Kenji Matsuda【序】π共役系分子ワイヤは、将来の分子エレクトロニクスにおいて導線の役割が期待されている分子の一つである 1。分子ワイヤの特性は、電子のトンネリングによる減衰定数 β で特徴づけられる。分子 1 つを流れるトンネル電流は長さ l に対して、量子化コンダクタンスを G0、減衰定数を β とすれば、G = G0 exp(−βl )の関係式で記述される。この挙動はトンネル電流のみに留まらず、
電子移動速度定数や交換相互作用にも見られることが実験により確認されている 2,3。これら電子移動や交換相互作用は共にπ共役鎖へのトンネル電流への寄与が期待されることから、量子化コ
ンダクタンスと同様の減衰定数を有する指数減衰が期待される。このことは電子移動もしくは交換相互作用の減衰定数を見積もることが出来れば、コンダクタンスの減衰定数を見積もることが出来ることを意味する。本研究ではこのような観点から、両端に有機ラジカルを結合させた分子について、ラジカル間の交換相互作用のπ共役系に対する長さに対する依存性を理論計算により検討した。
【計算】 ビラジカル間の交換相互作用の大きさは Heisenberg effective Hamiltonian
で記述される。交換相互作用の大きさは交換積分 J の大きさに比例する。山口らのグループはこの模型より交換積分 J を計算する理論式を導出した 4。

ここで、

,

はそれぞれ三重項、Broken-Symmetry (BS)状態のエネルギーであり、およびはスピン 2 乗演算子の平均値である。これは、三重項状態と一重項状態のエネルギー差を評
価し J を算出する方法である。これを BS 法と呼ぶことにする。計算に使用するプログラムは
Gaussian09 を、計算手法は DFT(UB3LYP/6-31g*および 6-311g**) を使用した。また、一重項・三重項状態を計算する際には基底状態である三重項状態の最適化構造を用いた。

【結果と考察】 我々はニトロニルニトロキシド(NN)とベルダジル(VER)を有機ラジカル種として選び、オリゴフェニレン(PP)とオリゴフェニレンエチニレン(PE)を π 共役ワイヤに選んだ。そして、以下の図にあるように NN-PP、NN-PE,VER-PP,VER-PE の 4 種類の分子群を設計し(図 1)、それらについて交換相互作用の距離依存性について調べた。

図 1．設計した分子群
それぞれの分子群に対して繰り返しユニットを 1 つずつ伸長させ、そのたびに交換相互作用を
ラジカル間の距離に対してプロットした。(図 2)3-1log (|2J| / cm )

4

2
1
0
-1
-2
0

10

20

30

40

Molecular Length / Å
図 2．交換相互作用の距離依存性○：VER-PP●：NN-PP△：VER-PE ▲：NN-PE

片対数グラフの傾きから、NN の系では PP に対して＝0.42 Å−1、PE に関しては＝0.24 Å−1 であることが計算された。コンダクタンスの実験では PP に対して β＝0.41Å−1, PE に対して β＝0.20Å−1 であることが報告されており 5,6、今回計算した値はこれをよく再現している。また、VERの系でも、PP に関して＝0.44 Å−1、PE に関して＝0.24 Å−1 であることが計算され、はスピン源である両端のラジカル種に依存せず、分子ワイヤを構成するユニット分子にのみ依存することが示された。これらの結果から、BS 法は有機分子ワイヤの値を予測する上で有用であることが実証された。

【参考文献】(1) A. Nitzan, M. A. Ratner, Science. 2003, 300, 1384. (2)K. Higashiguchi, K. Yumoto, K. Matsuda, Org.

Lett. 2010, 12, 5286. (3) B. Schlicke, P. Belser, L. De Cola, E. Sabbioni, V. Balzani, J. Am. Chem. Soc.
1999, 121, 4207. (4) K., Yamaguchi, Y. Takahara, T. Fueno, K. Nasu, Jpn. J. Appl. Phys. 1987, 26,
L1362.(5) Wold, D. J.; Haag, R.; Rampi, M. A. ; Frisbie, C. D. J. Phys. Chem. B, 2002, 106, 2813. (6) Lu,
Q.; Liu, K.;Zhang, H.; Du, Z.; Wang, X.; Wang F. ACS Nano 2009, 3,3861.

4P115歯車状両親媒性分子の自己組織化における理論的研究（横浜市立大院・生命ナノ 1, 東大院・総合 2, 産総研 3）
増子貴子 1, 山田健太 1, 平岡秀一 2, 長嶋雲兵 3, 立川仁典 1

Theoretical analysis of self-assembled gear-shaped amphiphile molecules(Yokohama city univ.1, Univ. of Tokyo2, AIST Institute3)
Takako Mashiko1, Kenta Yamada1, Shuichi Hiraoka2, Umpei Nagashima3, Masanori Tachikawa1【序】
親水部と疎水部をもつ両親媒性分子は、疎水部が水中で自己組織化することにより、ミセルや脂質二重層などの多様な会合体を形成することが知られている。この両親媒性分子およびその会
合体の会合の駆動力もしくは会合機構の詳細を解明することは、生物組織などの働きを理解する上で必要不可欠である。
近年、平岡らは図 1 に示すヘキサフェニルベンゼンに3 つの 3-ピリジル基を有する歯車状両親媒性分子が、水中にて一義的に六量体を形成することを見出した[1,2]。
ここで、置換基 R がメチル基の場合は会合し六量体を形成するが、他の置換基、例えば水素原子および水酸基のときは六量体を形成しない。このように、置換基 R からの異なる影響により、会合挙動が変化することが実験的に報告されている[1]。小関らは分子軌道法および分子動
力学（MD）シミュレーションにより、気相における歯車状両親媒性分子のメチル化体と水素化体の理論計算を実行した[3,4]。その結果、分子軌道計算による相互作
用エネルギーの観点から、π-π相互作用よりも CH-π相
互作用が分子間に与える影響が強く、メチル化体の方がより会合しやすいことが分かった[3]。また、会合体の温図 1 歯車状両親媒性分子の構造。置換基 R にメチル基、水素原子、水酸基がついたものをそれぞれメチル化
体、水素化体、水酸基体と表記する。

度の影響に注目し、温度を上げていくことによりメチル化体および水素化体のどちらも内部体積が大きくなり、さらにメチル化体よりも水素化体の方がその六量体構造を低い温度で保てなくなることが MD シミュレーションの結果から理解されている[4]。これらの研究により、水素化体よりもメチル化体が会合しやすいという実験の結果を支持する結果は得られているが、具体的に会合の駆動力に影響を与えている要因に関しては十分理解されていない。
そこで本研究では、置換基 R のみの違いで両親媒性分子の会合挙動が制御できるという点に着目し、MD シミュレーションを用い、歯車状両親媒性分子による六量体への自己組織化の要因を相互作用や構造変化の観点から解析する。

【計算方法】
AMBER9 プ ロ グ ラ ム [5]を 用 い 、 力 場 に は Generalized Amber Force F ield(GAFF) 、 電
荷 は HF/6-31G(d)に よ る 静 電 ポ テ ン シ ャ ル か ら 求 め ら れ た R ESP 電 荷 を 使 用 し 、 MD
シミュレーションを気相において実行した。また、六量体全体の相互作用を考慮
す る た め に 、 そ の カ ッ ト オ フ は 30.0 Å と し た 。

【結果・考察】
計算対象は図１に示した歯車状両親媒性分子の単量体および六量体であり、本研究では置換基 R としてメチル基、水素原子、水酸基の 3 種類を用いた。
実験で得られた X 線構造解析を参照構
造とし、まずは温度の揺らぎの会合体への 影 響 を 検 討 す る た め 、 NVT 計 算 を 行 った。そのメチル化体の六量体構造の最小化計算の結果を図 2 に示す。単量体および 六 量 体 の 分 子 構 造 に 関 し て 歯車状両親媒性分子の 1H-NMR の結果より、単量体では C3 対称、六量体ではその１つの構成分子が一種類であり、全てが C1 対称であるという示唆[2]がある。そこで、対称性を考慮しながら単量体の構造探索を行い、その構造・電荷を用いて六量体の新たな参照

図2構造を生成し、その六量体構造から NVT 計算を実車状両親媒性分子のメチル化体の六量体構

行し、置換基の違いによって現れる挙動の変化を造から最小化計算した構造。灰色が炭素、白議論した。特にメチル化体のみが会合体を形成すが水素、青が窒素原子である。

X 線構造解析の実験により得られた歯ることに注目し、MD シミュレーションの結果の解析から、歯車状両親媒性分子の六量体への会合、つまり自己組織化のメカニズムを決める要因
について検討を行なっている。
これらの詳細な結果に関しては当日発表を行う。

【参考文献】[1] S. Hiraoka, K. Hirano, M. Shiro, M Shionoya, J. Am. Chem. Soc., 130, 14368 (2008). [2] S. Hiraoka, K.
Hirano, T. Nakamura, M. Shiro, M. Shionoya, Angew. Chem., Int. Ed., 48, 7006 (2009). [3] J. Koseki, Y. Kita, S.
Hiraoka, U. Nagashima, M. Tachikawa, Theor. Chem. Acc., 130, 1055 (2011). [4] J. Koseki, Y. Kita, S. Hiraoka,
U. Nagashima, M. Tachikawa, Int. J. Quantum Chem., 10, 1002 (2012). [5] D. A. Case, T. A. Darden, T. E.
Cheatham, Ⅲ, et al., AMBER9, 9th ed., University of California: San Francisco, (2006).

4P-116制約密度汎関数理論を用いた PPV の三重項電子移動過程の解析（筑波大院・化）相川 小春，鞆津 典夫，守橋 健二

Constrained DFT Calculation on Triplet Electron Transfer Process of PPV(Univ. of Tsukuba) AIKAWA Koharu, TOMOTSU Norio, MORIHASHI Kenji【序】	
 Cao らは有機 EL 材料の一つである poly(para-phenylenevinylene) (図 1)の誘導体のelectroluminescence(EL) と photoluminescence(PL) 量子収率の比 QE(EL)/QE(PL) が 50%となることを示した1．これに対し，Shuai らは反応断面積の観点から説明を行った2．我々はこの現象を電子移動過程の速度定数の観点から説明することを目標とし，PPV の三重項電子移動過程の解析を

図 1: PPV

3

Constrained DFT を用いて行った．
【計算対象】	
 
PPV は ポ リ マ ー な の で モ デ ル 系 と し てtrans-stilbene (図 2), Oligo-PV の三量体
（図 3; 以下 OPV3），及び OPV3 の誘導体と

図 2: trans-stilbene

して両端にメチル基を置換したもの（図 4;
以下 dimethyl-OPV3）を用いた．
【計算手法】	
 上記の分子をそれぞれ，基底状態，三重項最図 3: OPV3低励起状態，アニオンラジカル，カチオンラ
ジカルの電子状態で B3LYP/6-31G(d)で構造最適化を行った．
この構造を以下のような組み合わせでダイ

図 4: dimethyl-OPV3

マーとして配置し，計算を行った．一つは始状態の電荷分離状態として D•−-A•+，もう一つは終状態のエキシトン状態として 1D-3A*の組み合わせである（D, A は計算上では同じ分子）．電子状態についてもそれぞれの構造において始状態，終状態の電子状態を再現するため CDFT を用いそれぞれ電荷，スピンに制約を 2.0 だけ適用した．
また，速度定数については以上のダイマー計算の結果を用いて Marcus 理論に基づく四点法で求めた．
【ダイマーの分子配向】	
 三つのモデル分子それぞれを平行に分子平面に対して 4 Å 離して配置した．また，

trans-stilbene についてはそれ以外の配向についても比較するため，ドナー分子について慣性主軸に沿って傾け・ずれをつくり計算を行った．そのほか，dimethyl-OPV3 については分子の結晶構造から傾きを持つ配向を抜き出し，そこにモノマー構造をはめたダイマーについて計算を行った．
【結果】	
 
平行配置については以下の表 1 のような結果となった．
表 1: 平行 配置の 各パラ メータ

ΔG˚ (kJ/mol)系

λ (kJ/mol)

Hab (kJ/mol)ket (/s)trans-stilbene−228.2

31.59−19.333.078×10−39

OPV3−162.725.15−19.401.561×10−18

dimethyl-OPV3−399.6253.4−22.381.125×1012

trans-stilbene, OPV3 では速度定数が非常に小さいが，dimethyl-OPV3 はそれらに比べかなり大きくなっている．また，trans-stilbene, OPV3 の二つで比較すると OPV3 のほうが 1020 ほども大きくなった．
また，慣性主軸に沿って傾け・ずれをつけた trans-stilbene はすべての結果で平行配置以下の速度定数となった．その中での傾向では，傾きのみの場合 30˚, 45˚, 60˚で比較すると 45˚が最もΔG˚とλの差が小さくなり速度定数が大きくなるという結果となった．
また，45˚の傾斜をつけた上で分子の伸長方向にずれをつけた場合，ΔG˚はずれが大きくなるほど小さくなるがλはほとんど変化しなかった．また 45˚の傾斜をつけた上で分子間距離に当たる垂直方向のずれを作った場合，近づけるほどΔG˚は大きくなりλはほとんど変化しなかった．
結晶構造にはめて計算した dimethyl-OPV3 の結果は以下の表 2 のようになった．T1,
T2 は結晶構造中の斜め配向のパターンである．
表 2: 結晶 構造の 配向で の各パ ラメー タ

ΔG˚ (kJ/mol)系

λ (kJ/mol)

Hab (kJ/mol)ket (/s)

T1−81.74

24.05

25.03

1.968×109

T2−331.426.08−11.201.225×10−141

Parallel−399.6253.4−22.381.125×1012

T1 ダイマーと前述の平行配置ではともに速度定数がでているが， T1 ダイマーではΔG˚が小さく，平行配置ではλが大きいために速度定数が大きくなっている．
1

Y. Cao et al., Nature 397, 414 (1999)
Z. Shuai et al., Phys. Rev. Lett. 84, 131 (2000)
3
Q. Wu and T. Van Voorhis, Phys. Rev. A 72, 024502 (2005)
2

4P117
分子動力学シミュレーションによるアルカリ金属ハライド水溶液の１次元振動分光スペクトル解析(京大院・理)伊藤広伸、谷村吉隆

Vibrational Spectroscopies of Aqueous Alkali Halide Solutions
in MD Simulation(Kyoto Univ. Science)

Hironobu Ito, Yoshitaka Tanimura【序】 水溶液中におけるイオンの水和のしくみを明らかにすることは、化学や生物における様々な現象を理解するうえで欠かすことができない。振動分光スペクトル解析は、
水和イオンの動力学的な振舞いを調べるのに適した方法であり、実験分野、理論分野を問
わず効果を発揮する。S. Meech らのグループは、様々な金属ハライド水溶液の低振動領域における等方性ラマン分光スペクトル解析を行い、濃度に依存して水と陰イオン間の水素結合モードの振動数が増加していくことを観測した [1]。この実験結果に興味を持ち、今回我々は水分子、イオン原子ともに分極ポテンシャルモデルを用いたアルカリ金属ハライド水溶液の分子動力学シミュレーションを行い、そのシミュレーション結果から１次元
赤外分光及びラマン分光スペクトルを計算することにより、水とイオン間の振動モードの解析を行った。
【計算方法、結果】 本研究では、水のポテンシャルモデルに当研究室で開発された分
極水モデル [2] を用いた。図１は、この分極水モデルの分子動力学シミュレーションにより計算した 300 K における純粋な水の１次元赤外振動分光スペクトルである。この図から分かるように、この分極水モデルは振動スペクトルにおいて低振動から高振動領域までの各振動モードを的確に捉えており、振動分光解析に適している。また、水とイオン間、
イオンとイオン間の相互作用に関したポテンシャルモデルを静電ポテンシャル、分散・反発ポテンシャル、分極ポテンシャルに分けて考えた。塩化ナトリウムをはじめとした強電解質水溶液中において静電相互作用はイオン雰囲気による静電遮蔽効果を受け、クーロ
ン力の及ぼす領域が非常に狭くなる。そのため、我々は静電ポテンシャルにデバイ遮蔽効
果を取り入れた湯川型ポテンシャル [3] を採用した。この静電相互作用について、水 − イオン間では単極子 − 単極子、双極子 − 単極子、四重極子 − 単極子、単極子 − 誘起双極子、双極子 − 誘起双極子、四重極子 − 誘起双極子、イオン − イオン間では単極子 − 単極子、単極子 − 誘起双極子、誘起双極子 − 誘起双極子までの静電相互作用を考慮した。また、分散・反発ポテンシャルには R. Vuilleumier らが第一原理計算により提案したポテ
ンシャルモデル [4] を用いた。イオン原子の分極については、用いた分極水モデルと同様に、誘起双極子モーメンに対して分極のエネルギーが最小になるような平衡条件 [5] により導出した。分極モデルを用いた分子動力学計算において、分子間距離が小さくなると、Polarization Catastrophe と呼ばれる数値的発散が生じることが知られている。これを防ぐため分子間距離が小さい領域において、静電ポテンシャルを補正する Thole 型 [6] の減衰関数を用いた。
本研究の分子動力学計算において、単極子 − 単極子、単極子 − 双極子 (誘起双極子)、双極子 (誘起双極子)− 双極子 (誘起双極子) の静電ポテンシャルは tinfoil boundary condition
によるエワルド法を用いて計算を行った。このエワルド和には Thole 型の減衰関数による補正を行っていない通常の湯川型ポテンシャルを用いた。なぜなら、減衰関数は分子間距離が小さい領域 (高々∼ 3 Å 程度) でしか影響を及ぼさないためである。時間刻みを 0.1
fs とし、運動方程式の数値的解法には、発展演算子の水の分子内振動に関する項に６次のシンプレティック差分法、その他の相互作用に関する項に２次のシンプレティック差分法を用いた。また、温度制御には能勢 − Hoover chain 熱浴を用いた。
実際の分子動力学計算の手順として、まず初めに共役勾配法による系全体のエネルギー最小化計算 (1000 step) を行い、NVT アンサンブルによる 5 K の極低温条件下での 100
ps のシミュレーションを行った。次に実験環境の 298.15 K まで温度を上げ、200 ps のシミュレーションにより平衡化を行った。１次元赤外分光及びラマン分光スペクトルの計
算には NVE アンサンブルによる 3 ns のシミュレーション結果を用いた。分子動力学計
算の解析の結果は当日発表する。

図 1: 300 K における水の赤外振動分光スペクトル【参考文献】[1] I. A. Heisler, K. Mazur, S. R. Meech, J. Phys. Chem. B 115, 1863 (2011).
[2] T. Hasegawa, Y. Tanimura, J. Phys. Chem. B 115, 5545 (2011).
[3] P. Debye and E. Hückel, Phys. Z. 24, 185 (1923).
[4] R. Vuilleumier et al, J. Chem. Phys. 136, 114507 (2012).
[5] J. Applequist, J. R. Carl, K.-K. Fung, J. Am. Chem. Soc. 94, 2952 (1972).
[6] B. T. Thole, Chem. Phys. 59, 341 (1981).

4P118光電子角度分布への配置間相互作用の影響（京大院理 1，CREST JST2，理研 3） 鈴木 喜一 1, 2, 3，鈴木 俊法 1, 2, 3Electron correlation effects in photoelectron angular distributions（Kyoto Univ.1，CREST JST2，RIKEN3） Y.-I. Suzuki1, 2, 3，T. Suzuki1, 2, 3電子相関は、物理と化学の様々な現象の中で重要な役割を果たしている。しかし、電子相
関が引き起こす現象を明確に研究する手段は限られている。これまで、二重イオン化で発生する二個の電子を「同時に」観測する方法などが試みられてきた。1 しかし、分子の対称性によっては、一電子放出過程でも研究が可能である。このことを明らかにするため、本研
究ではピラジン（点群 D2h）およびシクロプロパン（D3h）について電子相関を考慮した光電子
角度分布の計算を行なった。
ピラジンの S2→S1 非断熱遷移は、円錐交差を介した内部転換のベンチマークシス
テムである。この反応は、時間分解光電子画像観測法により、実時間で観測されている（時定数は 22±3 fs）。2, 3 実験では、光電子の運動エネルギー(PKE)ごとに、光電子角度分布を pump(264 nm)-probe(198 nm)の遅延時間の関数として測定した。光電子角度分布の変化は、電子状態が変化していることの証拠である。特に、PKE = 0.9 eV では S2→S1 の反応に伴う角度分布の変化がはっきりと観測された（図１a と 1b）。この PKE では、エネルギー的にイオンの電子基底状態(D0)のみが関与している。しかし、D0(n-1)←S2(ππ*) のイオン化過程は、独立電子近似では禁制であり、電子相関を考慮しなければ説明できない（図 2）。
図１cと 1d に電子相関を考慮した計算結果を示す。4 計算結果は実験を良く再現し、遅延時
間13 fsで観測された D0←S2 のほぼ等方的な光電子角度分布が良く説明できた。また、得られた波動関数を精査した所、電子相関は S2 と D0 のどちらの状態でも重要であることがわか
った。
13 fs

(b) > 100 fs

laser pol.

TRPEI obs.
(a)

(d)

D0 - S1

laser pol.

FOCI/CMSXα calc.
(c) D0 - S2

図 1: PKE=0.9 eV における、観測(a,b)および計算(c,d)された光電子角度分布

図 2: ピラジン S2 および D0 の電子配置

シクロプロパンは、160 nm 付近に強い吸収をもち、3p Rydberg (21E’) 状態に帰属されている。近接する一重項状態で、1E”, 1A1’, 1A2’ (3p Rydberg)状態は双極子遷移禁制である。3s Rydberg (11E’)状態は 180～190 nm 付近の弱い吸収に対応すると、考えられている。
独立電子近似では、3p Rydberg (21E’) 状態の Rydberg 軌道は e’、 3s Rydberg (11E’) 状態の Rydberg 軌道は a1’ 対称性を持っており、それらの軌道の対称性が光電子角度分布に反映する。5 一方、その二つの電子配置は同じ対称性 1E’ であり、電子相関を考慮した場
合の固有状態の波動関数は、それらの一次結合である。そのため、シクロプロパンの 1E’状態の光電子角度分布は、電子相関の影響を受け、a1’(3s)および e’(3p)軌道、両方の性質を反映すると予想できる。ここでは、3p Rydberg (21E’)状態に選択的に励起した場合の、共鳴二光子イオン化について計算を行なった。
図 3 に、計算で仮定したイオン化のスキームと、分子固定系での光電子角度分布を示した。4 励起、イオン化光の偏光は分子の三回軸（Z 軸）に、それぞれ垂直および平行とした。図 3c と 3d をくらべると、電子相関の影響が明らかである。電子相関を無視した場合の光
電子角度分布は、三回軸対称性を保持している（励起光の偏光方向には依存しない）。一方、電子相関を考慮することで、三回軸対称性が崩れている。この計算では、イオンは単一電子配置で良く記述されていて、対称性の崩壊は、電子相関による 3e’→3s Ryd.(14 %) と3e’→ 3p Ryd.(86 %)電子配置の混合が主な原因である。

図 3： (a)シクロプロパン分子についての、共鳴二光子イオン化のエネルギーダイアグラム。(b)空間に固定された分子に対する光の進行方向（青矢印）と偏光方向（両端矢印）。Θ は分子の二回軸と光の進行方向のなす角。(c)単一電子配置(3e’→3px,y)で計算された、分子固定系での光電子角度分布。(d)電子相関を考慮した場合の光電子角度分布。Θ = 110°, PKE = 0.1 eV。
1

例えば、T. Weber, H. Giessen, M. Weckenbrock, G. Urbasch, A. Staudte, L. Spielberger, O. Jagutzki,

V. Mergel, M. Vollmer, and R. Dorner, Nature 405, 658 (2000).
2

T. Horio, T. Fuji, Y.-I. Suzuki, and T. Suzuki, J. Am. Chem. Soc. 131, 10392 (2009).

3

Y.-I. Suzuki, T. Fuji, T. Horio, and T. Suzuki, J. Chem. Phys. 132, 174302 (2010).

4

計算方法： 中性(Φ0)およびイオン(Φ+)の電子波動関数は、配置間相互作用法(ピラジンは FOCI, シ

クロプロパンは CASSCF)で計算した。連続状態の一電子軌道(ψ)は、多重散乱法(CMSXα)で計算した。終状態の波動関数は、ψとΦ+の積を反対称化したもので定義した。分子構造は電子基底状態の平衡構造とした。
5

例えば、R. J. Smith, J. Anderson, and G. J. Lapeyre, Phys. Rev. Lett. 37, 1081 (1976).

緑色蛍光タンパク質変異体
CFP と GdFP に お け る 光 特 性 に 関 す る 理 論 的 研 究(横 市 大 院 生 命 ナ ノ 1 , 北 里 大 薬 2 ) 高 橋 麻 里 奈1 , 小 関 準 2,1 , 北 幸 海 1 , 川 島 雪 生 1 , 立 川 仁 典 1

Theoretical study on the absorption/fluorescence properties forthe mutants of green fluorescent protein, CFP and GdFP
(Yokohama City Univ.1 , Kitasato Univ.2 ) Marina Takahashi1 , Jun Koseki2,1 , Yukiumi Kita1 ,Yukio Kawashima1 , Masanori Tachikawa1【 序 】緑 色 蛍 光 タ ン パ ク 質 (GFP: Green
Arg96(+)
Fluorescent Protein) は 、1962 年 下 村 ら [1] に
Gln94
NH
O
よ り オ ワ ン ク ラ ゲ (Aequorea Victoria) か ら
H N
NH
HN
発 見 さ れ 、1990 年 代 に は そ の 遺 伝 子 の 単
R
Cro66
O
離 方 法 が 確 立 さ れ た 。そ れ 以 来 、細 胞 生 物
W6
His148
N
学や医学等の幅広い分野で生物マーカー
NN
Nと し て 利 用 さ れ て き た 。さ ら に 、現 在 に
1
H
NH
W
W1
O
W2
W5
至 る ま で 様々な 光 特 性 を も つ 蛍 光 タ ン パ
W3
HO
O
W4
ク 質 が 合 成 、発 見 さ れ て き た 。蛍 光 タ ン
OH
Thr203
パ ク 質 の 光 特 性 は 、主 に タ ン パ ク 質 の 内OHOO部 に 存 在 す る 発 光 素 (CRO: Chromophore) のSer205
Glu222(-)
構 造 に 依 存 す る こ と が 報 告 さ れ て い る [2]。
例 え ば 、GFP の CRO を 形 成 す る Tyr を Trp
Figure 1: The largest model for the active sitesに 変 異 さ せ た CFP(Cyan Fluorescent Protein)of CFP(R: H), (4-Me-Trp)CFP(R: CH3 ) andGdFP(R: NH2 ).
は 青 緑 色 蛍 光 を 放 ち 、さ ら に 、CFP の Trpに ア ミ ノ 基 を 付 加 さ せ た GdFP(Gold Fluorescent Protein)[3] は 、強 い 黄 色 蛍 光 を 放 つ 。
ま た 、多 く の 蛍 光 タ ン パ ク 質 に お い て 、CRO は 周 囲 残 基 や 水 分 子 と 水 素 結 合 ネット
ワ ー ク を 形 成 し て い る (Figure 1)。こ の 水 素 結 合 ネット ワ ー ク 変 化 は 、CRO の イ オ ン
化 状 態 変 化 を 誘 導 し 、光 特 性 に 変 化 を も た ら す こ と が 報 告 さ れ て い る [4]。
CFP の 吸 収 ス ペ ク ト ル は 、434, 452 nm に 弱 い ダ ブ ル ピ ー ク を も つ 形 状 と な る が 、
GdFP は 446 nm に シ ン グ ル ピ ー ク を も つ 。一 方 で 、GdFP の ア ミ ノ 基 を メ チ ル 基 置 換
し た (4-Me-Trp)CFP で は 、CFP と ほ ぼ 同 一 の ス ペ ク ト ル で あ る こ と が 報 告 さ れ て い
る [3]。さ ら に 、CFP に つ い て は 2 種 の 異 な る 立 体 配 座 (PDB ID: 1OXD, 1OXE) を 、ま た
GdFP に つ い て は 単 一 の 立 体 配 座 (1OXF) の X 線 結 晶 構 造 を 報 告 し 、異 な る ス ペ ク ト ル
形 状 は こ れ ら の 立 体 配 座 の 違 い が 原 因 で あ る と 説 明 し た [3]。こ れ に 対 し て Demachy
ら は 1OXD お よ び 1OXE に つ い て 分 子 動 力 学 計 算 を 行 い 、異 な る 立 体 配 座 に お け る 吸
収 エ ネ ル ギ ー に 違 い が な い こ と を 示 し た 。彼 ら は CFP に お け る 弱 い ダ ブ ル ピ ー ク はππ ∗ 及 び nπ ∗ 励 起 に よ り 異 な る 電 子 状 態 が 混 在 す る た め で あ る と 示 唆 し た [5]。
こ の よ う に 、CFP の 類 似 体 間 で 特 徴 的 な 光 特 性 の 違 い が 報 告 さ れ て い る [3] が 、その 光 特 性 変 化 の 原 因 に つ い て は 、未 だ に そ の 議 論 に 決 着 が つ い て い な い 。そ こ で 本
研 究 で は 、ア ミ ノ 基 置 換 に よって CFP か ら GdFP へ と 光 特 性 変 化 が 生 じ る 原 因 の 解 明
を 目 的 と し 、理 論 的 解 析 を 行った 。
2

2

【手法】

H

R

蛍 光 タ ン パ ク 質 の 光 特 性 は 主 に CRO 構 造に 依 存 す る こ と か ら 、π 共 役 に 関 与 す る 部

CRO66
H

H

終 端 さ せ る こ と で CFP と (4-Me-Trp)CFP お
る Model I)。構 造 は B3LYP レ ベ ル で 最 適 化
し 、励 起 状 態 計 算 は TD-B3LYP を 用 い て 1
点 計 算 を 行った 。基 底 関 数 は 全 て 6-31+G(d)を用いた。
ま た 、CRO 周 囲 環 境 を 考 慮 す る た め タ
ン パ ク 質 全 体 に 対 し て ONIOM(QM:MM) 計

H H

Val68
N

N

分 の 構 造 を 抽 出 し 、そ の 末 端 を メ チ ル 基 でよ び GdFP を モ デ ル 化 し た (Figure 2 に お け

O

H

N
H

[





N
H

T ru n c a te d p o s itio n ]: M o d e l I ( s u b s t . f o r -CH 3): M o d e l II (
-H )



HO

H

O
H

Leu64

N
CH3

H



H
O



Figure 2: The chromophore(CRO) sturucturesof CFP (R: H), (4-Me-Trp)CFP (R: CH3 ), and
GdFP (R: NH2 ). Model I and II are truncatedby methyl groups and hydrogens respectively.

算 を 実 行 し た 。QM 計 算 に は Model I と 同 一 の 計 算 レ ベ ル を 、MM 計 算 に は AMBER を
用 い 、PDB ID: 1OXD, 1OXF を 初 期 構 造 に 最 適 化 し た 。ONIOM 計 算 に お け る Model 領
域 は CRO 部 分 と 隣 接 す る 残 基 と の ペ プ チ ド 結 合 ま で 含 み 、水 素 で 終 端 さ せ た (Figure

2 に お け る Model II)。さ ら に タ ン パ ク 質 効 果 を み る た め に 、ONIOM 計 算 に よ る 最 適
化 構 造 か ら タ ン パ ク 質 部 分 を 除 去 す る こ と で 得 た Model II 分 子 に つ い て も Model I と同 様 に 励 起 状 態 計 算 を 行った 。

【結果と考察】Table 1 に 計 算 に よって 得 ら れ た 吸 収 エ ネ ル ギ ー と 、対 応 す る 実 験 値 を 示 し た 。CFPと (4-Me-Trp)CFP の 吸 収 ス ペ ク ト ル に 関 す る レッド シ フ ト の 実 験 値 は 0.02-0.03 eV で ある が 、気 相 Model I に よ る 計 算 結 果 は 0.07 eV と な り、実 験 値 を 定 性 的 に 再 現 し て いる 。同 様 に 気 相 モ デ ル に よ り、GdFP の レッド シ フ ト も 定 性 的 に 再 現 で き た 。そ の 他の 手 法 に よ る 計 算 結 果 も CFP と 比 較 し た (4-Me-Trp)CFP お よ び GdFP の レッド シ フ トに 関 し て 定 性 的 に 実 験 値 と 一 致 し た 結 果 を 得 る こ と が で き た 。ま た 、Model II を 用
い た 気 相 及 び タ ン パ ク 質 中 に お け る 計 算 結 果 を 比 較 す る こ と に よ り、タ ン パ ク 質 効
果 は CRO の 吸 収 ス ペ ク ト ル を 0.04-0.07 eV ブ ル ー シ フ ト さ せ る こ と が わ かった 。
さ ら に 、CFP と GdFP に お い て ス ペ ク ト ル 形 状 の 違 い が 生 じ る 原 因 を 明 ら か に す
る た め に 、ONIOM 計 算 に よ る 最 適 化 構 造 を 用 い て CRO 周 囲 に あ る 各 残 基 あ る い は
水 が 励 起 エ ネ ル ギ ー に 及 ぼ す 影 響 に つ い て 解 析 を 行って い る 。こ れ ら の 結 果 詳 細 については当日ポスターにて発表する。
Table 1: Absorption energies for CFP, (4-Me-Trp)CFP and GdFP. [eV]
CFP
(4-Me-Trp)CFP GdFP
Gas phase
In Protein
Exptl.

Model I
Model II
Model II

3.31
3.10
3.14

3.24
3.06
3.13

2.73
2.84
2.90

2.74, 2.86

2.72, 2.83

2.66

[1] O. Simomura, F. H. Johnson and Y. Saiga, J. Cell. Comp. Physiol., 59, 223 (1962) [2] R. Y. Tsien, Annu.
Rev. Biochem., 67, 509 (1998) [3] J. H. Bae, M. Rubini, G. Jung, G. Wiegand, M. H.J. Seifert, M. K. Azim, J.
Kim, A. Zumbusch, T. A. Holak, L. Moroder, R. Huber and N. Budisa, J. Mol. Biol., 328, 1071 (2003) [4] Noam
Agmon, Biophys. J., 88, 2452 (2005) [5] I. Demachy, J. Ridard, H. L. Pasquier, E. Durnerin, G. Vallverdu, P.
Archirel and B Lévy, J. Phys. Chem. B, 109, 24121 (2005)

4P120多原子分子の全自由度を考慮した解離性電子付着断面積の計算方法の開発（埼玉大院理工） 杉岡雄仁 , 高柳敏幸

Development of a computational method of
dissociative electron attachment cross sections forpolyatomic molecules in full dimensionality(Saitama Univ.) Yuji Sugioka , Toshiyuki Takayanagi
低エネルギー電子が分子に衝突すると、一時的にラジカルアニオンが生成し、結合が解離する。
この過程は解離性電子付着と呼ばれ、プラズマ化学や放射線化学等の様々な分野で極めて重要な役割を果たしている。解離性電子付着の起こりやすさは電子の衝突エネルギーに対する解離性電子付着断面積で表されるが、実験的に絶対値を決めることは極めて難しく、理論的に算出する方法が開発されてきた。しかし、核の全自由度を量子論的に扱うことは分子の大きさが増えるほど非常に困難であるため、これまでの理論は一部の自由度のみを考慮していた[1]。しかし、これでは解離による生成物のエネルギー分配は分からない。よって、核の全自由度を考慮した計算方法の開発が必要である。
本研究では経路積分分子動力学法と共鳴複素ポテンシャル理論を組み合わせた計算方法を開発した。初期構造と初期運動量を求める際には、中性分子のポテンシャルエネルギー曲面上で経路
積分分子動力学法を用いて、温度と核の量子性を含めたサンプリングを行なった。次に、得られた初期構造と初期運動量をもとにアニオン分子のポテンシャルエネルギー曲面上で古典トラジェ
クトリー計算を行い、アニオンの寿命を積算し、解離性電子付着断面積を求めた。この方法を実験・理論ともに研究されている H2O (H2O + e → OH + H ) [2,3]と CF3Cl (CF3Cl + e → CF3 +Cl ) [1,4]に適用し、計算方法の妥当性を検証した。解離性電子付着断面積は温度に大きく依存する系があるため[4]、その温度依存性を再現できるかどうかについても考慮した。
Fig.1,2 はそれぞれ H2O と CF3Cl における解離性電子付着断面積の電子の衝突エネルギー依存性を表している。H2O については厳密な量子動力学計算で求められた値とよく一致し、CF3Cl については実験値を定性的に再現している。また、Fig.3,4 はそれぞれ H2O と CF3Cl の解離における H,CF3 の相対運動エネルギー分布の電子の衝突エネルギー依存性を表している。核の全自由度を扱った計算をすることによって、このように解離性電子付着における詳細な情報を得ることができる。

Fig.1 電子の衝突エネルギーに対する H2O

Fig.2 電 子 の 衝 突 エ ネ ル ギ ー に 対 す るの解離性電子付着断面積

CF3Cl の解離性電子付着断面積

Fig.3 電子の衝突エネルギーに対する H

Fig.4 電 子 の 衝 突 エ ネ ル ギ ー に 対 す るの相対運動エネルギー分布

CF3 の相対運動エネルギー分布[1] M. Tarana, K. Houfek, J. Horá ček, I. I. Fabricant, Phys. Rev. A, 84 (2011) 052712.
[2] D. S. Belić , M. Landau, R. I. Hall, J. Phys. B 14 (1981) 175.
[3] D. J. Haxton, T. N. Rescigno, C. W. McCurdy, Phys. Rev. A 75 (2007) 012711.
[4] I. Hahndorf, E. Illenberger, L. Lehr, J. Manz, Chem. Phys. Lett. 231 (1994) 460.

4P121
4P121
統一的なフェルミ準位を用いない分割統治型量子化学計算法の開発（早大先進理工1，早大理工研2，JST-CREST3）

吉川 武司1，中井 浩巳1, 2, 3

Development of quantum-chemistry calculation based on the divide-and-conquer method
without common Fermi level(Advanced Science and Engineering, Waseda Univ.1, RISE, Waseda Univ.2, CREST, JST Agency3)
Takeshi Yoshikawa1, Hiromi Nakai1, 2, 3
【緒言】
緒言】
これまでに大規模電子状態計算を可能とするために、フラグメント分子軌道(FMO)法[1]やエロンゲーション法[2]等の様々な分割型計算理論が開発されてきた。なかでも分割統治(DC)法[3,4]は部分系に対して統一的なフェルミ準位を用いることで非局在系や開殻系を初めとした
様々な分子に適用が可能であり、この統一的なフェルミ準位が DC 法の利点の一つであった。
しかし、原子数が増加するにつれ、統一的なフェルミ準位を決定するときに多くの時間を要するため、SCF 計算の並列化において、大きなボトルネックとなっている。一方で、分子間相互作用で集積している分子系では、各分子における電子数はほとんど変化しないため、フェルミ準位を決定することなく、高精度かつ高効率に計算することが可能である。そこで、本発表では、統一的なフェルミ準位を決定しない方法（DC-nofermi 法）の結果を示し、分子性結晶、溶媒系や分子間相互作用系等の分子集合系への有用性を提示する。
【方法論】
方法論】
DC法は全系をいくつかの部分系に分けて計算を行うことで計算コストを削減する方法である。部分系の周りのバッファ領域を含めた局在化領域で部分系αの軌道を構築することで、周囲の効果を取り込むことができる。
DC-SCF 計算においては、以下の制約条件(1)を課し統一的なフェルミ準位εF を決定することにより、各局在化領域内の電子数を自動的かつ適切に振り分けることができる。
αne = ∑ pµν∑ fβ (ε F − ε iα )CµaiCναi Sνµ
α

(1)

i

ここで p は分割行列、µ, νは原子軌道、i は分子軌道、fβはフェルミ関数、εαは部分系の軌道エ
ネルギー、Cαは部分系の軌道係数、S は重なり行列を示している。フェルミ関数を用いることで、非整数状態の電子数も取り扱うことが可能となるため、電子が非局在化した系に対しても精度よく計算することができる。
DC-nofermi 法では、各局在化領域の電子数をあらかじめ決定して、部分系の計算を行う。
水やフッ化水素等の分子集合系では各分子間での電子移動が小さく、電子は分子内に局在化しているため、統一的なフェルミ準位を用いなくても高精度に計算することが可能である。
この方法を用いる利点としては、統一的なフェルミ準位を決定しないため、部分系の Loop
が 2 回から 1 回に削減されることによって通信効率を高め、通常の方法よりも並列化効率が高くなることが挙げられる。

SCF

・全体の F 行列の生成・全系の F 行列から局在化領域の Fα行列を抜き出す
DC-SCF
Loop subsystem α・Fα行列の対角化
End Loop
・統一的なフェルミ準位準位の決定
Loop subsystem α・電子数と軌道係数を用いて Dα行列を求める
End Loop

DC-nofermi
Loop subsystem α・Fα行列の対角化・電子数と軌道係数を用いて Dα行列を求める
End Loop・分割行列を用いて全体の D 行列を求める
・エネルギー、電子密度等のプロパティを求める

Fig. 1. Scheme of DC-SCF and DC-nofermi method

【数値検証】
数値検証】
新たに開発した DC-nofermi 法を用いて、フッ化水素(FH)20 量体の計算を行った。DC 法による計算では、部分系は FH 1 分子からなるユニットとし、バッファは左右 4 ユニットとした。
Table 1 には、様々な基底関数に対する全エネルギーと、従来法と DC 法とのエネルギー誤差を示した。DC 法と DC-nofermi 法のどちらの場合でも、6-311++G**を除けばそのエネルギー
誤差は 0.3 mhartree 以下と非常に精度よく計算することに成功した。また、通常 DC 法において分散関数を用いる場合、統一的なフェルミ準位付近に軌道が多く存在するため、フェルミ
準位を決めることが困難となる。その結果、適切にフェルミ準位を決定することができず誤差が大きくなる傾向がある。一方で、DC-nofermi 法では統一的なフェルミ準位を決定する必要がないため、通常の DC 法より高精度に計算することが可能となる。

Basis-set
6-31G
6-31G**
6-31++G**
6-311G
6-311G**
6-311++G**

Table 1. Basis-set dependence of the total energies for (FH)20(diff.)
DC-HF(diff.)
DC-nofermi-HF[hartree][mhartree][hartree][mhartree]-1999.974958(+0.001)-1999.974958(+0.001)-2000.487735(+0.000)-2000.487735(+0.000)-2000.667559(+0.069)-2000.667427(+0.201)-2000.715097(+0.167)
Not conv.
-2001.152997(+0.115)-2001.153126(-0.015)-2001.241676(-1.567)-2001.242648(-2.540)

Conv.
[hartree]-1999.974959-2000.487735-2000.667627-2000.715264-2001.153111-2001.240108

次に、He216、(H2)1000、(FH)20 と(H2O)125 に対する全エネルギーと計算時間について Tabel 2と Table 3 に示す。基底関数は 6-31G とする。DC 法による計算では、部分系は 1 原子/分子とし、バッファは 5Åとした。どの分子集合系においてもそのエネルギー誤差は 1.5 mhartree 以
下と小さく従来法を精度良く再現することに成功した。また、DC-nofermi 法を用いることによって従来の DC 法よりもさらに高速に計算することに成功した。このことから、DC-nofermi法を用いることによって、分子集合系については、従来の DC の精度を大きく損なうことな
くかつ高速に計算が可能であることがわかった。当日は、巨大分子が入り混じった分子集合系の取り扱いについても報告する予定である。

He216(H2)1000(FH)20(H2O)125

Table 2. Total energies of He216, (H2)1000, (FH)20 and (H2O)125
DC-nofermi-HF(diff.)
DC-HF(diff.)[hartree][mhartree][hartree][mhartree]-616.714914(0.000)-616.714914(0.000)-1125.423726(-0.241)-1125.423652(-0.315)-1999.974958(+0.001)-1999.974958(+0.001)-9364.889363(-1.442)-9364.889432(-1.311)Conv.
[hartree]-616.714914-1125.423967-1999.974959-9364.888131

Table 3. CPU time and SCF iteration of He216, (H2)1000, (FH)20 and (H2O)125
DC-nofermi-HF
DC-HF
Conv.
CPU time
CPU time
CPU time
iteration
iteration
iteration
[min]
[min]
[min]
He216
0.6
3
0.6
3
0.6
3
(H2)1000
112.9
4
128.0
4
149.4
4
(FH)20
0.1
9
0.2
10
0.2
10
26.5
11
44.5
18
288.5
23
(H2O)125[1] K. Kitaura, E. Ikeo, T. Asada, T. Nakano, and M. Uebayashi, Chem. Phys. Lett. 313, 701 (1999).
[2] A. Imamura, Y. Aoki, and K. Maekawa, J. Chem. Phys. 95, 5419 (1001).
[3] W. Yang, Phys. Rev. Lett. 66, 1438 (1991).
[4] M. Kobayashi and H. Nakai, in Linear-Scaling Techniques in Computational Chemistry andPhysics: Methods and Applications (2011, Springer), pp. 97-127.

4P122残差切除法による Coupled Perturbed 方程式の解法（東大生研 1，ASI 総研 2）○阿部 敏彦 1，関根 義人 2，佐藤 文俊 1
【はじめに】
楕円型偏微分方程式の数値解法として優れた収束性を有す残差切除法を、Coupled Perturbed (CP)方程式の解法に適用した。本方法は、CP 方程式の解法においても優れた収束性を示し、小さい分子については Pople 法と同等な性能を示した。分子のサイズが大きくなるにつれて、収束到達までのステップ数が増えるが、Pople 法と同等か、もしくは勝る結果を与えた。分子サイズが大きな系では、残差切除法の特徴から、特に収束性において本方法が有利となる可能性が示された。

【方法】
GIAO 法に基づく原子軌道に関する CPHF 方程式は、占有軌道と空軌道への写像表現として、(1)と表される[1]。ここで とた、はそれぞれ Hartree-Fock (HF)方程式を満たす LCAO 係数と電子密度行列である。まはそれぞれ外部磁場に関する 1 次の摂動密度行列、摂動 Fock 行列および摂動重なり行列である。CPHF 方程式を解いて
程式は

を得る。

は

に依存しており、その係数をすべて展開すると、CPHF 方を未知数とする線形方程式となり、(2)と表せる。ここで はの要素を並べてベクトルとして表現したものである。 は行列要素数が軌道数の 4 乗となるスーパーマトリクスであり、通常式(2)は繰り返し計算で解かれる。
残差切除法による計算の手順は次のようになる[2]。
(1) 初期値 を与える。
(2) 残差を求める。
(3)を、内部ソルバによる最小単位の反復で暫定近似解を求める。
(4) 合成摂動量による近似解の残差最小化係数
を
計算する。
(5)
から新しい近似解を得る。
(2)から(5)の手順を繰り返すことにより、残差を減少させていく。この残差切除法を CP 方程式へ適用するにあたり、本研究では(3)において暫定近似解を求める代わりに、
(3)
とした。

【実験結果および考察】
残差切除法と Pople 法[3]について、Hartree-Fock 法による GIAO 法に基づく NMR 遮蔽テンソル計算で用いる
CPHF 方程式の繰り返し解法[4]の収束計算実験を行った。残差切除法ではとの 2 種類の値の結果を示した。比較のため、単純な繰り返し計算で解いた simple damping 法も一部実行した。テスト用分子として水分子、
窒素分子、ペプチド鎖を用いた。ペプチド鎖は高度好熱性細菌の putative acylphosphatase の全体の 87 残基から、32～37 残基目を 6 残基ペプチド鎖、1～28 残基目を 28 残基ペプチド鎖としてそれぞれ切り出した構造を参考に初期構造を作成した。基底関数は 28 残基ペプチド鎖には 6-31G を、それ以外の分子には 6-31G*を用いた。以下に 6
残基ペプチド鎖と 28 残基ペプチド鎖での結果を示す。
残差切除法の収束性能が Pople 法を顕著に上回っている。これは、Pople 法では、繰り返しベクトル集合の強い
線形従属性のため、直交化の精度が保てずに、収束が停滞したためである。また、6 残基よりも分子サイズが大きい 28 残基の場合は、この問題はより顕著に表れており、そのため到達可能な最小の残差ノルムは一桁以上大きくなっている。

図 1. 6 残基ペプチド鎖計算における、Pople 法、残差切除法(L=5, L=10)、simple damping 法の各手法による残差ノルムの収束過程。繰り返しステップ数(横軸)に対する残差ノルムの 10 の対数(縦軸)を示す。右図は分子構造。

図 2.

28 残基ペプチド鎖計算における、Pople 法、残差切除法(L=5, L=10)の各手法による残差ノルムの収束過程。

キャプションは図 1 と同様である。右図は分子構造。

本研究は「イノベーション基盤シミュレーションソフトウェアの研究開発」により行われた。計算には、東京大学情報基盤センターの計算機を使わせて頂いた。

参考文献[1] M. A. Freitag, B. Hillman and A. Agrawal, J. Chem. Phys., 120, 1197, 2004.
[2] A. Tamura, K. Kikuchi and T. Takahashi, J. Comp. Phys., 137, 247, 1997.
[3] Pople, Raghavachari, Schlegel and Binkley, Int. J. Quantum Chem. Symp., 13, 225, 1979
[4] R. Ditchfield, Mol. Phys., 27, 789, 1974.

4P123
スレーター行列式を使用した量子モンテカルロ法：
F12 法との組み合わせによる高精度計算（神戸大システム情報） ○大塚 勇起, 天能 精一郎

High-accuracy calculations by a hybrid F12 and determinantal-based
quantum Monte Carlo method（Kobe Univ.） Yuhki Ohtsuka and Seiichiro Ten-no【序】 拡散モンテカルロ法（プロジェクタモンテカルロ(PMC)法）は、高精度かつ並列化効率の高い理論である。スーパーコンピューターでの計算に適しているため、近年、注目を集めている。この方法では、電子は粒子として扱われる。波動関数の反対称性を考慮するためには、あらかじめ計算された試行波動関数（ガイド関数）の節が利用されるが、その結果として、計算精度は使用されたガイド関数に強く依存するという問題が生じる。
我々は、PMC 法のガイド関数依存性を回避するために、電子を粒子ではなくスレータ
ー行列式(SD)用いて表す方法を提案した(PMC-SD 法)[1,2]。スレーター行列式によって波動関数の反対称性が自動的に満たされるので、ガイド関数は不必要となる。この方法では、
Full-CI 解を近似できるため、結合解離など既存の方法では再現が困難な系にも応用が可能であると考えている。しかしながら、拡散モンテカルロ法では、極限として数値的厳密解が得られるのに対して、PMC-SD 法の精度は、スレーター行列式を生成するために用いられた基底関数に依存してしまう。近年、基底関数の遅い収束性を改善するために、電子間の距
離を波動関数に含める方法（F12 法もしくは、R12 法と呼ばれる[3,4]）が急速に発展しているが、本研究の目的は、PMC-SD 法と F12 法を組合わせることにより、適度な大きさの基底関数を使用してシュレーディンガー方程式の精密な解を求めることである(PMC-SD-F12 法と呼ぶ)。
【理論】 PMC-SD 法では、虚時間 τ における波動関数は、スレーター行列式（電子配置）
I と、 I であるウォーカーの個数 nI (τ ) を用いることによって以下のように定義される。
Ψ PMC-SD (=τ)

1∑ ( nI (τ ) )2

× ∑ nI (τ ) =II∑ C (τ ) I
I

II

また、スレーター行列間の遷移確率や電子エネルギーは以下のように計算できる。
電子配置 I から J への遷移確率

I (1 − ∆τ H ) J

虚時間 τ における電子エネルギー

EPMC-SD (τ ) =ここで、 0 には、通常 Hartree-Fock 配置を選ぶ。

0 H Ψ PMC-SD (τ )C0 (τ )

PMC-SD-F12 法では、遷移確率や電子エネルギーを計算するときに、通常のハミルトニアン
H の代わりに、F12 法によって相似変換されたハミルトニアン=
H ′ exp(− R) H exp( R) を使用する。ここで、R は以下のように定義される。
1
R = ∑ αβ Q12 R12 ij Eα i Eβ j
2 α ,βEα i は、unitary group generator であり、 R12 と f12 は、それぞれ論理推進演算子とスレーター型

ジェミナル相関因子を表す。
31=
αβ R12 ij
αβ f12 ij + αβ f12 ji ,
8
8
1
f12 =
− exp ( −γ r12 )γ

今回は、強い直交化演算子 Q12 として Ansatz 1 と呼ばれるものと、より精度が高くスタンダー
ドな Modified Ansatz 2 の 2 種類を使用した。

( P =ϕ (1) ϕ (1) )Q12 =−(1 P1 )(1 − P2 )Ansatz 1

Modified Ansatz 2 Q12 =(1 − O1 )(1 − O2 ) − V1V2

1

p

p

( O =ϕ (1) ϕ (1) , V =ϕ (1) ϕ (1) )
1

i

i

1

a

a

【結果と考察】 TABLE 1 に、PMC-SD 法と PMC-SD-F12 法を Ne 原子に応用した結果を示す。PMC-SD-F12 法によって計算された相関エネルギーは、PMC-SD 法の結果よりも基底関数に対して速く収束した。Modified Ansatz 2 を使用した結果は、Ansatz １の結果と比較して、さらに収束性が良い。今回の最も精度の高い計算（Modified Ansatz 2, aug-cc-pCVQZ）
では、現在、最高精度の結果を与える CCSDT-R12 法[5]（390.6(4) mEh）と、ほぼ同じ結果(390.6(9) mEh)が得られた。発表当日は、アルゴリズムと計算の詳細、他の系への応用例を紹介する予定である。
PMC-SD 法と PMC-SD-F12 法による Ne 原子の相関エネルギー (in mEh)

TABLE 1

aug-cc-pCVXZ
MethodCCSD(T)
PMC-SD
PMC-SD-F12 (Ansatz1)
PMC-SD-F12 (Ansatz2)

D

T

Q

256.5
256.3(6)315.0(4)381.6(7)344.9344.6 (8)365.3 (9)388 (1)372.0
372 (1)382.0(9)390.6(9)Y. Ohtsuka and S. Nagase, Chem. Phys. Lett., 463, 431, (2008).
Y. Ohtsuka and S. Nagase, Chem. Phys. Lett., 485, 367, (2010).
S. Ten-no, Chem. Phys. Lett., 398, 56, (2004).
W. Klopper, F. R. Manby, S. Ten-no, and E. F. Valeev, Int. Rev. Phys.
Chem. , 25 427-468 (2006).
[5] T. Shiozaki, M. Kamiya, S. Hirata, and E. F. Valeev, J. Chem. Phys., 130,054101-1-10 (2009).

[1]
[2]
[3]
[4]

４P-124
緑色蛍光タンパク質の光励起動力学：QM/MM 法による非断熱 ab initio 分子動力学(上智大学 1,京都大学福井謙一記念研究センター2) 北淑宏 1, 南部伸孝 1, 石田俊正 2
Photoexcitation dynamics of Green Fluorescent Protein : nonadiabatic ab initio MD
approach with QM/MM scheme(Sophia University1, Fukui Institute for Fundamental Chemistry, Kyoto University 2)
Yoshihiro Kita1, Shinkoh Nanbu1, Toshimasa Ishida2【序】緑色蛍光タンパク質(Green fluorescence protein : GFP)はオワンクラゲ内でイクオ
リンと複合体を形成し、p-ヒドロキシルベンジリデンイミダゾリジノン発色団（図 1）がβ-バレル内に囲まれている安定な分子である。単独でも励起光さえあれば発光することか
らマーカータンパク質としての活用が期待され、医療など様々な分野への利用が期待されている。しかしながら光励起された GFP の詳細な発光過程はまだ解明されていない。その一方、発色団のみを考慮した計O

算では無輻射遷移過程となる非断熱遷移を伴う経路や、発色団と周N

囲のアミノ酸におけるプロトン移N動（図 1）が考えられている。

O

以前の我々の研究では、発色団

非断熱遷移などを観察した。本研

R2

Hとその周辺のアミノ酸構造をもとに、励起状態動力学計算を実行し

H

O Wat25
H
H
O

究では発色団を含む全タンパク質

HO
O
O
R3

を考慮した分子動力学計算により

R4

Ser205

光吸収後の電子励起状態ダイナミ
クスを解明する。

R1

Glu222

図 1 発色団と周囲のアミノ酸とのプロトンワイヤー【理論計算】QM/MM(ONIOM)モ
デルを用い、非断熱 ab initio 分子動力学計算を行う。QM 層には 、発 光 に 関 わ る と 予
想 さ れ る 発 色 団 と 、 水 分 子 25、 セ リ ン 205 、 グ ル タ ミ ン 酸 222 を 選 ん だ 。
MM 層 に は 残 り の β-バレル構造を持つ分子を選んだ。量 子 化 学 計 算 の レ ベ ル は
CASSCF/MIDI と し た 。 な お 、 今 回 は 分 子 埋 め 込 み (mo lecular embedding :ME 法 )を 用 い た た め 、MM 部 の 電 荷 に つ い て は 考 慮 に 入 れ ず に 計 算 を 実 行 した 。MM 計算の部分は Amber の ff03 力場を用いた。Wigner 分布に従い、光吸収後の発色団の部分の初期構造と運動量を乱数で決定した。速度ベルレ式を用いて得られた初期構
造と運動量を用い、核座標を時間発展させ、特に各タイムステップにおいて上記量子化学計算を実施し、断熱ポテンシャルエネルギーから原子核の運動を求めた。さらに、ポテン
シャルエネルギーの擬交差が生じた場合には Zhu-Nakamura 理論を採用した古典軌道ホップ（TSH）法により非断熱遷移を考慮した。

【結果と考察】励起状態動力学計算を行った２つのトラジェクトリーで異なるプロトン移動が確認された。
1 つ目のトラジェクトリーでは光励起後 588 fs まで計算を進めた。ポテンシャルエネル
ギーに関しては、S2 と S1 の間で何度も擬交差が確認されたが、非断熱遷移は起きずに S1状態にとどまっていた。しかし S2 と S1 の断熱ポテンシャルが接近しているため、S2－S1間の非断熱遷移の可能性は示唆される。また、プロトン移動に伴う O-H 核間距離については、光吸収後 94 fs に協奏的なプロトン移動が見られた。ここで、プロトンの移動には
O-H 結合距離（0.96Å）を参考にした。しかしその後、3 つのプロトンは元の位置に戻り、その付近で振動を繰り返した。
2 つ目のトラジェクトリーでは、255 fs まで計算を進めた。ポテンシャルエネルギーに関しては上述のトラジェクトリーと同様の挙動を示した。また、プロトン移動に伴う O-H
核間距離の変化では、光励起後 70 fs で協奏的なプロトン移動が見られた。その後元の位置に戻らず、その付近で結合状態を保ちながら振動を繰り返した。一般的にプロトン移動
はピコ秒（1000 fs）単位で起こるとされていたが、今回はその 10 分の 1 の時間で確認することが出来た。
発色団と周囲のアミノ酸のみを考慮した、以前の我々の研究では 100 fs までに非断熱遷移が多く確認されていたが、本研究では現時点でまだ観測されなかった。理由としては周りのタンパク質の存在が発色団の振動を妨げているため、非断熱遷移時に生じるような五員環と六員環の cis-trans 異性化に起因するような分子振動が起こりにくいことが考えられる。今後は、量子化学計算手法の改善と、発色団以外の部分の電荷を考慮した electron
embedding(EE)での計算を行うことで、励起状態プロトン移動や電子配置を正確に描写した計算を行っていく予定である。

図 2 観測された励起状態プロトン移動

4P125

化学ポテンシャルによる多電子ダイナミクス解析:
LiH 分子の強レーザー場によるイオン化（東北大院理 1, 横市大院生命ナノ 2, 東大院理 3, 大阪府大院理 4）

大村周 1, 小山田隆行 2,

河野裕彦 1, 加藤毅 3, 小関史朗 4

A Chemical Potential Analysis of Multielectron Dynamics:
Intense-Field Induced Ionization of the LiH Molecule(Tohoku Univ.1, Yokohama City Univ.2, The Univ. of Tokyo, Osaka Prefecture Univ.) Shu Ohmura1,
Takayuki Oyamada2, Hirohiko Kono1, Tsuyoshi Kato3, Shiro Koseki4
【序】 近年のレーザー技術の発達により、高出力(光強度～1020 W/cm2)かつ超短(フェムトからアト秒)のレーザーパルスが得られるようになってきた。このようなレーザー光を分子に照射すると、トンネルイオン化や高次高調波発生といった興味深い現象が観測される。最近では、近赤外強レーザー場において HOMO よりも下の軌道からのトンネルイオン化を示す実験
報告もなされ注目を集めている[1]。このような非摂動的な電子のダイナミクスを求めるための第一原理計算手法として、我々は多配置時間依存 Hartree-Fock (MCTDHF) 法を開発してきた[2][3]。さらに詳細な解析のために、非線形、非定常の問題にも適用できる時間依存の化学
ポテンシャル𝜇(𝑡)などを導入している[4]。H2 分子の場合、各軌道からのイオン化効率とその軌道の𝜇(𝑡)の上昇が相関しており、電子ダイナミクスを理解する上で時間依存の化学ポテン
シャルが有用であることが明らかになってきた。本研究では化学ポテンシャルによる解析の応用として内殻を持つ最小系である LiH 分子の強レーザー場に対する応答を計算し、𝜇(𝑡)を用いた解析を行った。
【理論】 MCTDHF 法では、多電子波動関数Φ(𝑡)は時間依存分子軌道で構成されるスレータ
ー行列式Φ𝐼 (𝑡)の線形結合
Φ(𝑡) = ∑ 𝐶𝐼 (𝑡)Φ𝐼 (𝑡)(1)𝐼で表される。ここで𝐶𝐼 (𝑡)は時間依存の展開係数(CI 係数)である。この波動関数に Dirac-Frenkelの時間依存変分原理を適用することによって、軌道関数と CI 係数の両に対する運動方程式が導出される。係数だけでなく軌道の最適化も行うことで、より少ない配置数で電子波束の時
間発展を追うことができる。また開発した MCTDHF 法では数値軌道関数を用いるので、束縛状態と連続状態を等価に扱うことができ、イオン化のような大振幅の現象を記述することが可能となっている。この MCTDHF 法を用いて、𝜇(𝑡)を計算した。自然軌道 j の化学ポテンシ
ャル𝜇𝑗 (𝑡)は、その軌道の占有数 wj(t)と全エネルギー𝐸(t)によって𝜇(𝑡) ≡ ∂𝐸⁄𝜕𝑤𝑗 と定義され、𝐸(𝑡) = ∑ 𝑤𝑗 (𝑡)𝜇𝑗 (𝑡)(2)𝑗という関係が成り立つ。この𝜇𝑗 (𝑡)を用いて自然軌道ごとのエネルギーの時間変化を定量化することができる。

【結論と考察】 本研究では分子軸 z 偏
光の 2 サイクル電場 (波長 1520 nm、最大強度 I ～3.0×1012 W/cm2)と相互作
用する LiH 分子をとりあげる。1、2、
3軌道を用い、9 配置考慮した。Li→H

電場の方向を電場の正の向きとしている。図1 は各時刻での系のエネルギー変化と双極子変化を示している。電場が切れた時
点で全エネルギーは上昇しており、励起

図 1：トータルエネルギーの変化量(赤線)と状態の混入がみられる。エネルギーや双誘起双極子変化(青線)、及び用いた電場極子の変化は初めの約 1.5 サイクルは断(破線、arb. unit)

熱的だが、徐々に電場に対して非断熱的な応答を示すようになる。このようなふるまいをより詳しくみるために(3)式に従って全エネ
ルギーを分解し、軌道ごとの寄与を調べた。その結果，各軌道の双極子モーメントが初めの
約 1.5 周期の間，電場に対し 断熱的に変化していることがわかった。この様子は各軌道の時
間発展のスナップ ショットでも確認できる。 たとえば，2 軌道のスナップショットをみると (図 2(a)-(c) )、分子軌道が電場の変化に応じて左右に振られていることがわかる。このような双極子の変化に対して、軌道の化学ポテンシャルの変化はよい相関を示す。すなわち化学ポテンシャルが、電子相関を伴う電子の動きと対応して変化しており、多電子ダイナミク

10

(a)

電場からの力の向き

(b)
H

Li

5
0

bohr

15

スを理解する重要なツールとなることがわかった。詳細は当日発表する。

-20

-10

0
z (bohr)

10

20
(c)

図 2：分子軸 z 偏光のレーザー電場と相互作用する LiH 分子の自然軌道()の絶対値の常用対数プロット。円筒座標を用いている。それぞれ時刻(a)t =0 fs、(b)t =4.5 fs、(c)t =7.0 fs での様子。

-6【参照論文】[1] H. Akagi et al., Science 325 1364 (2009).
[2] T. Kato and H. Kono, Chem. Phys. Lett. 392 533 (2004).
[3] T. Kato and H. Kono, J. Chem. Phys. 128, 184102 (2008).

[4] T. Kato and H. Kono, Chem. Phys. 366, 46 (2009).

-4

-2

0

2

4P126
実空間密度汎関数法による表面吸着系の光学応答量の計算手法の実装(阪大院基礎工) 乾 智也，奥野 克樹，馬場 剛史，米田 京平，福田 幸太郎，
重田 育照，中野 雅由

Theoretical study on optical properties of molecules adsorbed on surface
based on real-space density functional theory(Osaka Univ. Graduate School of Engineering Science) Tomoya Inui, Katsuki Okuno,
Takeshi Baba, Kyohei Yoneda, Kotaro Fukuda, Yasuteru Shigeta, Masayoshi Nakano【序】分子にレーザー光のような強い光が照射されると、分子に分極が誘起されるが、その分極は光電場のべき級数に展開できる。この高次の項に起因して起きる光学現象を非線形光学応答といい、3 次元光メモリなどへ応用が期待されている。これまで我々の研究室では大きな第二超分極率を持つ、非線形光学応答特性の優れた様々な分子の設計を行ってきた。従来の研究ターゲットは単一分子や有限サイズの分子クラスターが主であったが、新たなターゲットとして環境の効果に着目している。ここでいう環境とは溶液内などの均一場の効果[1]ばかりでなく、表面やケージ内など不均一場の効果も含み、それらが非線形光学特性にどのような影響を与えるかは非常に興味深い。光学応答は電子のゆらぎである分極に起因しているため、たとえば分子がある表面系に吸着したような系では電荷移動が起こり光学応答特性の変化が誘起される。これにより、より優れた光学応答特性をもつ分子の設計、及び光学応
答特性の制御が期待される。
	  固体表面を対象とする際には二次元の無限系を取り扱う必要がある。表面平行方向に対しては周期境界条件を用いることで取り扱うことができるが、表面垂直方向には周期性は存在しない。このような系を取り扱う方法として、クーロンカットオフ法[2]を実空間密度汎関数（RSDFT）法に実装する。RSDFT 法では空間を有限のグリッドに区切り、有限差分を用いる
ことで Kohn-Sham 方程式を解くことから、並列化が容易であるといった利点をもつ。
	
  本研究では表面系での線形・非線形光学応答特性の評価のために、二次元の周期性を持つ系に対する RSDFT 法プログラムを実装し、固体表面に対して結果の妥当性を検証する。
【理論】クーロンポテンシャルの場合、逆空間でのポテンシャルは 3 次元周期境界条件下で以下のように表される。

V (G) = n (G) v 3D (G),!!

v 3D (G) =

4!
G2

(

)

ここで n ( G ) は電荷密度のフーリエ変換であり、 G = G x , G y , Gz は逆空間ベクトルである。これに対してクーロンカットオフ法では 2 次元周期系の場合に周期性の無い方向（ここでは z軸方向）のクーロン相互作用が単位胞内で打ち切られるように変更を加えることで 2 次元系を取り扱う。具体的に、上式の v 3D ( G ) は 2 次元周期において、次のように変更される。

{

' 4!
!Gxy ( Z /2 )" Gz / Gxy sin (Gz Z / 2 ) ! cos (Gz Z / 2 )$)) G 2 1 + e#%
2Dv (G) = (4!
) 2 1 ! cos (Gz Z / 2 ) ! (Gz Z / 2 ) sin (Gz Z / 2 )
G
)* z

(

)

{

}

(a)

} (G & 0 )xy(G = 0 )xy(b)

ここで、Gxy = Gx2 + Gy2 、Z は単位胞の z 軸方向長さである。RSDFT 法のプログラム[3]に上記の２次元クーロ
ンカットオフ法を実装した。また、２次元系の z 軸方向に対する分極率・超分極率の取り扱いとして、有限 y場法を用いた計算法を実装した。ここで、分極率を有限場法で求めるために、z 軸方向に外部電場の掛かっx
図 1. 単位胞の概形

た計算のためのプログラム変更も同時に行った。
【計算】本研究では交換相関汎関数に PZ81 を採用し、(a)グラフェン (Lx=4.648, Ly=8.050, Lz=20.0)

図 1 の単位胞（赤枠）を用い、xy 方向には通常の周期(b)GNR(Lx=4.648, Ly=16.10, Lz=20.0)

境界条件を、z 軸方向（紙面垂直方向）には前述の方法を用いることで、2 次元グラフェン（以後グラフェ
ンと呼ぶ）(a)、一次元グラフェン（グラフェンナノリ
ボン; GNR）の１種であるポリアセン(b)の(超)分極率を算出した。またグラフェン及び GNR との比較のためにそれぞれ C6H6、C24H12、C54H18 及びポリアセン系（ベンゼン環数 2−12）の分極率についても算出した。
こ の 計 算 は 交 換 相 関 汎 関 数 に PW91 、 基 底 関 数 は6-31G*を採用し GAUSSIAN09 を用いた。
【結果と考察】図 2 にグラフェン及びその類似分子の

図 2. グラフェンとその類似分子の1 炭素原子あたりの分極率

1 炭素原子あたりの分極率を示す。ベンゼンから順に炭素の数が増えていくに従って 1 炭素原子あたりの分極率はある値に収束していくが、本プログラムの結果はその値に近い値が得られていることがわかる。図 3
は GNR における結果を示している。GNR においても前述のグラフェンの場合と同様にポリアセン系の分極率によく一致した結果が得られている。また 1 炭素
原子あたりのグラフェン、ポリアセンの第二超分極率
γ はそれぞれ 22.1 および 81.6 a.u.であった。ただしグ

図 3. GNR とその類似分子のラフェンのγの計算では計算精度などの理由から、単位1 炭素原子あたりの分極率

胞を xy 方向に二倍に拡大した 16 原子の単位胞を用い ( n はポリアセン系におけるベンゼン環の数)た。表面吸着系に関する結果は当日報告する。
【参考文献】[1] 乾智也、他、第５回分子科学討論会 [4P105] (2011). [2] C. A. Rozzi et al. Phys. Rev. B 73,
205119 (2006). [3] J.-I. Iwata et al., J. Comp. Phys. 229, 2339 (2010).

4P127
分子内のラグランジュ面による化学結合の理論的研究(京都大学 1 )埜崎寛雄1 、市川和秀 1 、立花明知 1

Theoretical study of chemical bond by Lagrange surface(Kyoto Univ.1 )Hiroo Nozaki1 、Kazuhide Ichikawa1 、Akitomo Tachibana1

近年、電子ストレステンソル τ Skl [1](k, l = 1, 2, 3) (定義を下の式 1 に示す) と言う、物質内部において電子に加わる力と関連付けられるテンソルに対し、これを用いた化学結合に対する理論的研究が行われている。例として、τ Skl の発散により定義されるテンションと言うベクトル場に対し、それがゼロとなる点をラグランジュ点と言う結合を代表する点として考える研究 [2] などを挙げる事が出来る。これは、ラグランジュ点におけるエネルギー密度 (電子ストレステンソルのトレースにより得られる) を比較することで結合次数を計算するものであった。さらに、ラグランジュ点におけるストレステンソルの固有値の縮退および最大固有値の符号から、金属結合と共有結合を識別し得る可能性も、昨今では示唆されている。このように従来、電子ストレステンソル、およびそれから派生するテンションは、化学結合に対する理論的研究に対して用いられてきた。
本研究では、このテンションに対し、そのセパラトリクスを分子を原子毎の固有の領域に分割する境界面として考えた。このテンションのセパラトリクスが、ラグランジュ面 [3] と呼ばれる面である。このラグランジュ面によって、分子内部における原子毎の固有の領域が求まれば、領域内の電荷の積分によって、分子内で原子が持つ電荷を定義する事ができる。
最後に、具体的なラグランジュ面の形状を、HF、HCl、F2 と言った 3 つの二原子分子 (電子状態は全て singlet) の結合軸を通る面に対し、次頁の図 [1,2,3] に示す。ここで、図中で太線で示されている部分
こそがラグランジュ面であり、矢印はテンションベクトルの向きである。テンションのセパラトリクス、
即ちラグランジュ面によって、分子内における原子固有の領域が明らかにされている事が確認できる。
なお、各原子は図の白丸の位置に存在し、原子核の質量が大きい原子が左側に存在する。これら分子に対する構造最適化および計算は、プログラムパッケージ Gaussian09 を用い、計算方法 CCSD、基底関
数系 CC-pVQZ により行った。実線は分子の表面を表すと考えられている electronic interface(S) [3]
であり、これは電子の運動エネルギー密度 nT (⃗r) [3] が 0 となる領域である。
本発表では、このラグランジュ面によって原子固有の領域を求めその内部の電子密度を成分する事で得られる原子毎の電荷と、従来から存在する Mulliken 電荷との比較を行う事を予定している。
なお、先に述べた電子ストレステンソルと電子の運動エネルギー密度は、次に示す式 1、式 2 によって与えられる。

τ

Skl

[
]
~2 ∑
∂ 2 ψi (⃗r) ∂ψi∗ (⃗r) ∂ψi (⃗r) ∂ 2 ψi∗ (⃗r)∂ψi∗ (⃗r) ∂ψi (⃗r)∗(⃗r) =
νi ψi (⃗r)−+ψi (⃗r) −.(1)
4m i∂xk ∂xl∂xk∂xl∂xk ∂xl∂xl∂xk

nT (⃗r) = −~2 ∑νi [ψi∗ (⃗r)∆ψi (⃗r) + ∆ψi∗ (⃗r)ψi (⃗r)] .
4m i

(2)

これらの式において、~ は換算プランク定数であり、m は電子の質量である。また ψi は自然軌道、νi
はその占有数である。

HF

HCl
1

0.8

1
1

0.8
0.6

0

0.4

−0.4
−0.8
−0.8 −0.4

図1

0 0.4
x[Å]

0.4

0.2

−0.5

0.2

0

−1

0

0.8

−1 −0.5

1
1

0.8

0.5

0.6

0

0.4

y[Å]

0 0.5
x[Å]

図 2 HCl 分子のラグランジュ面

F2

−0.5

0.2

−1

図3

0.6

0

HF 分子のラグランジュ面

0−1 −0.5 0 0.5
x[Å]

0.8

0.5
y[Å]

y[Å]

0.4

1

F2 分子のラグランジュ面

参考文献[1] A. Tachibana, J. Chem. Phys. 115, 8 (2001).
[2] P. Szarek, A. Tachibana, J Mol Model. 13, 651 (2007).
[3] A. Tachibana, J. Mol. Structure. 943, 138 (2010).

1

4P128対称性からみた原子クラスターのモード間エネルギー移動と反応機構(早大基幹理工)

岡 右里恵, 柳尾 朋洋

Symmetry Dependence of Vibrational Energy Transfer andthe Mechanisms for Structural Transitions of Atomic Clusters(School of Fundamental Science and Engineering, Waseda Univ.)【はじめに】

Yurie Oka, Tomohiro Yanao

原子分子集合系の自己組織化や集団運動の仕組みを理解することは，現代の分子科学における重要な課題である．一般に，原子分子集合系の集団運動は，非線形性，非平衡性が強く，また高次元であるため，従来の振動モード解析法を適用する上でいくつかの困難を有する．
例えば，構造変化を伴うような大振幅の集団運動においては，モード間の強い非線形結合のため，
従来の基準振動解析の適用は困難である．また，従来の基準振動解析では，Eckart フレーム[1]を用いて振動と回転の相互作用（Coriolis 相互作用）[2]を近似的に無視するが，実は，この種の相互作用は大振幅運動においては無視できない重要な役割を果たしている可能性がある．
そこで本研究では，超球座標の手法[3]をもとに，超球モード解析[4]と呼ぶ新たなモード解析法を発展させ応用する．超球モード解析の利点は，系の質量分布の変化に注目することで，系の構造変化を少数の集団変数の運動に縮約できる可能性がある点にある．また，超球モード解析では，
Eckart フレームではなく慣性主軸フレームを用いるため，系の変形と回転の相互作用の効果を近似せずに正しく扱うことができる．本発表では，ナノ秩序構造体のプロトタイプとしてしばしば用いられる 6~13 原子の Lennard-Jones クラスター（以下 Ar クラスターと呼ぶ）の構造転移運動を取り上げ，その力学的機構を探求する．特に，系の自発的な構造変化の引き金を引く分子内エネルギー移動過程と，構造変化の動的な駆動力の役割に焦点を合わせる．
【超球モード解析】 超球モード解析においては，n 体系がもつ 3n − 6 個の変形（振動）の自由度は，3 つの回転半径モード，
3 つの捩りモード，および 3n −12 個のずりモードによって記述される．例として，PBP と呼ばれる Ar! 7 クラスターの１つの局所平衡構造の近傍における上述の３種類のモードを図１€

図１: Ar7 クラスターの PBP 構造近傍における(a)回転半径モード a1 ，(b)捩りモード
ω 23 ，(c)ずりモード γ14 の例．

に示した．3 つの回転半径モードは，系の瞬間的な３本の慣性€主軸の方向への質量の伸縮運動を表し，各回転半径モードに対€€応する回転半径を大きい順に a ,a ,a と表す．また，3 つの捩
1

2

€

3

(ij = 12, 23, 31) とし， 3n −12 個のずりモードの角速度を γ ik (i = 1,2,3,k = 4,,n −1) と表す．図２に全エネ

りモードの角速度を ω ij

ルギー一定，全角運動量ゼロのAr7 クラスターが PBP から€€
COCT と呼ばれる構造へと変化する際の3 つの回転半径€€a1,a2 ,a3€の時間変化を示す．図２より，回転半径は系の質量分€

布の対称性の変化に応じて値を大きく変える集団変数の役割

図２: Ar7 クラスターの PBP 構造からを果たしていることが窺える．実際，構造変化がまさに発生しCOCT 構造への構造転移における３つのている時間帯（図２の２本の点線で挟まれた時間帯）において，回転半径 ai (i = 1,2,3) の時間変化．

クラスター内の全モードへの運動エネルギーの分布を解析す
ると， a2 モードと a3 モードが他の多くのモードに比べて特に大きな運動エネルギーをもって支配€€的な役割を果たしていることが分かった．これは，図２に見るように，PBP と COCT の間の構造変化においては２つの回転半径 a2 ， a3 が特に大きく変化することと対応している．

!

!
!

!

【構造転移の力学的機構】 では，回転半径 a2 ， a3 にエネルギーを注入し，系の構造変化を引き起こす駆動力は何であろうか？この問いに答えるには，回転半径に作用する力を解析することが有効である．一般に， n 体系の３つの回転半径 a1,a2 ,a3 に対する Euler-Lagrange 方程式は，
n−1
d 2 ai22!
! 2 ∂V ,
=
Γ
ω
+
Γ
ω
+
a
ij ij
ik ik
i ∑ γ im −
dt 2
∂ai
m= 4

(i, j,k ) = (1,2,3), (2,3,1), (3,1,2)!ω = −ω ）．ここで，右辺第４項は，原子間相互作用のポテンシャルに起因する通常の形となる（€ijji€のポテンシャル力である．一方，右辺第１項，第２項は捩りモード
ω ij ，ω ik に起因して回転半径に€作用する「捩りの遠心力」を表す． Γij は ai と a j に依存する係数である．一方，右辺第３項の和に€含まれる各項は，ずりモード
γ im に起因して回転半径に作用する「ずりの遠心力」である．これら€€２種類の遠心力は系の捩りおよびずり運動に起因する内部遠心力であり，角運動量に起因して生じる通常の遠心力とは異なる．

€€
図３には，Ar7 クラスターの構造変化の直前と直後の時間
1.0（無次元）の間に回転半径に作用する全内部遠心力と全ポテンシャル力の変化を，時間間隔 0.125 ごとに示した．
図３において，負の時間領域
は構造変化前，正の時間領域
は構造変化後を表す．また，

図３: Ar7 クラスターの構造転移の直前，最中，直後において回転半径 ai

時刻 0 と示された２本の垂直は力の強弱を表し，赤色は正の力，青色は負の力を表す． Fi ≡ Fdc,i + Fpot,i は全内線に挟まれた時間帯（上部にT と記載）は，構造変化の集(i = 1,2,3)に作用する全内部遠心力 Fdc,i および全ポテンシャル力 Fpot,i の時間変化．色の濃淡
部遠心力と全ポテンシャル力の合力を表す．(a)は PBP から COCT，(b)は COCT か€ €ら PBP へと構造転移を行った多数の軌道についてのデータを平均化している．

€

団運動がまさに発生している€€

時間帯を表す（図２の２本の点線で挟まれた時間帯に対応）．この図を詳しく分析すると，PBP
から COCT への構造変化の直前（図３(a)の時刻 t = −0.125 ）において， a2 の負方向に強い内部遠心力 Fdc,2 が作用し，a3 の正方向に強いポテンシャル力 Fpot,3 が作用していることが分かる．つまり，これらの力が駆動力となって，図２で見たように， a2 は小さく， a3 は大きくなる方向へと加速さ€れ，構造変化が開始するのである．一方，COCT から PBP への逆反応においては，構造変化が!

!

発生する直前（図３(b)の時刻
t = −0.125 ）において， a2 の正方向に強い内部遠心力 Fdc,2 が作用し，!

€a3 の負方向に強い内部遠心力 Fdc,3 が作用していることが分かる．つまり，これらの力が駆動力と!
!
なって， a2 は大きく， a3 は小さくなる方向へと加速され，構造変化が開始するのである．
€

!

【まとめ】

!

!

以上のように，本研究ではまず，原子クラスターの構造転移を実質的に支配する少!

!数の集団変数を，モード間のエネルギー移動の観点から抽出した．さらに，これら集団変数を適!
切な方向に加速して，系の構造変化を引き起こす駆動力の一端を明らかにした．ポテンシャルに由来する静的な力のみではなく，内部遠心力のような動的な力が系の構造変化の開始に大きく関与している点が重要である．本研究の手法を各種のナノ秩序構造体の構造転移運動に適用し，その駆動機構を明らかにすることは，今後の重要な目標である．

【参考文献】[1] C. Eckart, Phys. Rev. 47, 552 (1935). [2] X. Chapuisat and A. Nauts, Phys. Rev. A 44, 1328
(1991). [3] R. G. Littlejohn and M. Reinsch, Rev. Mod. Phys. 69, 213 (1997). [4] T. Yanao, W. S. Koon, and J.
E. Marsden, J. Chem. Phys. 130, 144111 (2009).

4P129ポリヒドロキシフラーレンの近赤外光による振動励起と反応ダイナミクスに関する理論的研究(東北大院理 1 名古屋大院理 2)池田旭伸 1，新津直幸 1，菅野学 1，Stephan Irle2，河野裕彦 1

Theoretical study on vibrational excitation and reaction dynamics
of polyhydroxy fullerene induced by near-infrared light(Tohoku Univ.1 Nagoya Univ.2) Akinobu Ikeda1, Naoyuki Niitsui1, Manabu Kanno1, Stephan Irle2, Hirohiko Kono1【序】水酸化フラーレンC60(OH)24に定常レーザー(785 nm)を照射すると、低強度(<103 W/cm2)でも高熱を持ち衝突や連続的な崩壊を経て、カーボン
ナノチューブや多層のフラーレンになることが知られている(Fig. 1)[1]。カーボンナノチューブは優れた強度や弾性をもち、構造により伝導性が変化するため、高強度材料や半導体などへの応用が期待されており、この反応はナノチューブの新しい合成法として興味深い。フラーレンC60や他の官能基を持つ誘導体では同様の反応は見られず、詳しい反応機構は解明されていない。この反応の初期段階を調べるため、定常レーザー中によって加熱された水酸化フラーレンC60(OH)24の挙動を理論
計算により評価した。

Fig. 1 水酸化フラーレン C60(OH)24 から
ナノチューブへの反応過程【手法】定常レーザー中の C60(OH)24 の挙動を調べるために DFTB(Density Functional based
Tight Binding)法によりダイナミクス計算を行った。DFTB 法とは、ヒュッケル法のような緊密結合法に密度汎関数理論を近似的に取り入れた高速で精度の良い電子状態計算を行う手法である。
用いたレーザーは 785 nm の波長を持つ矩形波である。原子核の運動は、各時刻のレーザー電場
を静電場として求められた時間依存断熱ポテンシャル面上の古典粒子として記述する時間依存断熱状態法[2]により評価した。低強度レーザーによる
長時間ダイナミクスの計算は困難であるため、本研
究 で は 光 強 度 を 109, 1011, 1013 W/cm2 と し てC60(OH)24 に積算される電場の効果を評価した。さらに、光エネルギーにより振動励起した様子を調べるために，C60(OH)24 とより簡単なモデル系である
C60(OH)4 に高い余剰エネルギーが注入されたと仮定して，その後のダイナミクスを計算した。 両方の計算で見られた脱水反応を調べるため、
B3LYP/6-31G**レベルの密度汎関数法(DFT)により
C60(OH)2 の遷移状態計算を行った。

Fig. 2 C60(OH)24 の OH 伸縮振動の周波数
スペクトル (I = 1013 W/cm2)

【結果】定常レーザー中の C60(OH)24 において水酸基の振動励起がみられた。伸縮振動の強度

Table 1 Fig. 2 に見られる周波数スペクトル現れる２本のピークの光強度依存性を見積もるために、水酸基の結合長の時間変化 I(W/cm2)

E (V/m)Peak 1

Peak 2(0-2 ps)をフーリエ変換し、スペクトルを計算し 1013

8.68×109

2.34×100.811×103

た(Fig. 2)。すると、伸縮振動に相当する~3800 1011
cm-1 のピーク (Peak 1)とレーザーの周109

8.68×108

2.17×1040.816×104

8.68×107

2.90×1050.815×105

3

波 数 に 相 当 す る 12740 cm-1 の ピ ー ク(Peak 2)が得られた。さらに、ピークの光
強度依存性を見てみると(Table 1)、Peak
2 の高さが電場の強さ E に比例している
ことから水酸基が非共鳴電場の力により強制振動していることが分かる。Peak 1
に関しても電場に依存して大きく変化す
ることが読み取れる。
さらに C60(OH)24 の運動エネルギーから求めた温度の時間変化を出力したとこ
ろ、~10ps の周期で大きな振幅が観測された(Fig. 3)。この振幅は分子全体の回転に対応している。C60 ではこのような反応は見られないため、分子全体の双極子モーメント

Fig. 3 定常レーザー中の C60(OH)24 の温度変化(I = 1013 W/cm2)が関係しているのではないかと予想される。当日の発表ではより詳しい解析の結果を発表する。
高熱状態の C60(OH)24(~2000 K)において隣り合った水酸基から脱水反応が起こることが確認された。この反応を詳しく調べるため、隣り合った水酸基を二箇所もつ C60(OH)4 のダイナミクスを計算したところ、同じ反応機構と考えられる脱水がみられた。まず片方の水酸基と炭素の結合長(RC-O)が伸び(~2.0Å)、
もう片方の水酸基の水素が隣の酸素へ移動し水分子を形成する(Fig. 4)。水分子の脱離後、水素を失った酸素は一時的に直Fig. 4 C60(OH)4 の脱水反応の様子

立し(~100 fs)、隣の単結合や二重結合に倒れこむことにより安定化する。倒れる方向により生成物がそれぞれエーテルとエポキシドとなり生成比は 5：3 程度であった。DFT 計算からポテンシ
ャル曲面を求めると、RC-O の伸長によりもう片方の水酸基の水素が移動しやすくなることが分かり、この機構は妥当といえる。C60(OH)4 のアレニウスプロットと DFT の遷移状態計算から脱水反応の活性化エネルギーは生成物によらず 2 eV 程度と見積もられた。C60 の解離の初期段階として起こる C2 脱離の活性化エネルギーは約 10 eV であることが知られており[3]、得られた値はこれに比べて十分に小さい。このことから、実験においても反応の初期段階として脱水が起きていることが予想される。
[1] V. Krishna, N. Stevens, B. Koopman, and B. Moudgil, Nature Nanotechnology 5, 330 (2010).
[2] Y. Sato, H. Kono, S.Koseki, and Y. Fujimura, J. Am. Chem. Soc. 125, 8020 (2003).
[3] S. Matt, O. Echt, P Scheier, and T. D. Märk, Chem. Phys. Lett. 345, 194(2001).

4P130Rigged QED に基づいた原子核場の量子統計性とその挙動の関係の理論的研究（京大院工）高田 崇二郎，宮里 敏秀，池田 裕治，瀬波 大土，立花 明知 1

Theoretical study of the relation between the statistics of nucleus fields and their behavior
based on Rigged QED(Department of Micro Engineering, Kyoto University) Sojiro Takada, Toshihide Miyazato, Yuji Ikeda,
Masato Senami, Akitomo Tachibana【研究の背景・目的】
近年、量子力学的現象が巨視的なスケールで観測できる Bose-Einstein 凝縮（BEC）が脚光を浴びている。BEC は電子数が同じ原子であっても原子核がボソンかフェルミオンかにより、BEC を起こすかどうかが決まる。
また、QED（quantum electrodynamics）に基づく研究は自由粒子に対しては完成度の高いものとなっているが、束縛状態に対しては満足のいくものとなっていない。
本研究では、光子の放出・吸収が相互作用の本質的役割を果たすという立場に立ち、原子核も場の演算子として取り入れた立花の提唱する Rigged QED[1],[2] に基づき、束縛状態系の時間発展を計算するコードを開発し、原子核の統計性による影響を調べる。
【理論】
Rigged QED の中でも電子を２成分場で表したものを Primary Rigged QED と呼ぶ。本研究では非相対論の範囲内での議論を行う。非常に重い原子や非常にエネルギー状態の高い状況を考えない
限りは非相対論の範囲内の議論で十分であるが、ラグランジアンの中に全ての展開項を用いればローレンツ対称性を回復することが知られている。また、電子・原子核について以下のようなハ
ミルトニアンを用いた。
電子：

原子核：
実際の計算には Heisenberg 描像に基づき計算を行っており、任意の時間における演算子を、初期状態における演算子で展開している。種々の演算子は時間発展に伴い、展開に用いた演算子の次数が爆発的に増大するため、高次の演算子に対しては期待値を取ることで対処した。
また、原子核がボソンかフェルミオンかでそれが満たす交換関係・反交換関係は異なり、種々の物理量においてその統計性による違いを調べる。
【結果及び考察】
初期状態の原子の中心を点(x,y,z) = (0,0,±2.55) [bohr]に置いた２つのリチウム原子に対して、計算を行った。図１、２に原子核がボソン、フェルミオンそれぞれの場合についての結果をまとめ、それぞれについて初期状態における原子核密度と時刻 0.024、0.06、0.12[a.u.]における原子核
密度の変化量を示した。
原子核がボソンの場合は図１(b)~(d)から、原子核の中心から外側へと密度が流れていく傾向が見て取れる。逆に、原子核がフェルミオンの場合は図２(b)~(d)から、原子核の中心へと密度が流れていく傾向が見てとれる。このように丁度反対方向への運動を引き起こす効果は交換相互作用
によって作りだされたものと考えられる。
【今後の計画】
本研究では Rigged QED に基づき原子核がボソンかフェルミオンかの違いによりその挙動の違いを確認する事ができた。今後は、原子核を展開する際に用いる基底関数として、どのような関数系を用いるべきか考察する。
また、現段階では交換関係・反交換関係が時間発展するにつれて、初期状態で満たしていた関係からずれを生じてしまっている。この事はどのように再規格化を行って行くべきかという議論と密接に関係している。
現在行っているのは全電子数・全原子核数の保存に基づいた再規格化であり、今後、物理量の再規格化をさらに改良する方法について模索する。
(a)

(b)

(c)

(d)

図１ 原子核がボソンの場合(a) 原子核密度の初期状態(b)~(d) 左から時刻 0.024, 0.06, 0.12[a.u.]における原子核密度の初期状態からの変化量
(a)

(b)

(c)

(d)

図２ 原子核がフェルミオンの場合(a) 原子核密度の初期状態(b)~(d) 左から時刻 0.024, 0.06, 0.12[a.u.]における原子核密度の初期状態からの変化量【文献】[1] A. Tachibana, J. Mol. Struct. : THEOCHEM 943, 138 (2010).
[2] A. Tachibana, J. Mol. Modeling 11, 301 (2005).
akitomo@scl.kyoto-u.ac.jp1