4P077
ドデカンチオールを用いたイオン液体相からヘキサン相への金ナノクラスター抽出における抽出時間効果（千葉大院融合 1、日大文理 2）
大西 慧 1，畠山 義清 2，加藤 淳一 1，西川 恵子 1
Extraction Time Effect in the Extraction of Gold Nanoclusters to Hexane Phasefrom the Ionic Liquid Phase using Dodecanethiol(Chiba Univ.1, Nihon Univ.2)
Kei Onishi1, Yoshikiyo Hatakeyama2, Jun-ichi Kato1, Keiko Nishikawa1【序】イオン液体(IL)は常温において液体状態で存在する有機塩であり、難揮発性や高いイオン伝導性といった性質を持つ物質である。特に難揮発性を利用し、スパッタ法により、IL 中に金ナノクラスター(AuNCs)を調製する手法が報告されている[1]。我々は本手法により得られる AuNCsの粒径について制御条件検討を行い、粒径制御が可能であることを明らかにしてきた[2-4]。スパッタ法は多くの利点を有する調製手法であるが、少なからず粒径に分布を持つことが課題である。
この課題を解決するため、水－トルエン系において AuNCs をサイズ選択的に抽出する報告[5]に着目した。
同様に IL 中に存在する AuNCs においても抽出によってサイズ選別ができないかと考え、種々の粒径分布を有する IL 中の AuNCs について、1-Dodecanethiol を用いたヘキサン相への抽出実験を行った。小角 X 線散乱測定、UV-Vis 吸収測定による情報から、AuNCs のサイズによる抽出時間依存性について報告を行う。
【実験】これまでの研究経験より、AuNCs 調製に用いる IL を、粘度が低いなど比較的扱いやすい

IL である 1-butyl-3-methylimidazolium tetrafluoroborate とした。AuNCs 調製はスパッタ法により行った。このとき試料温度を 20,40,60,80 °C と変化させることで粒径分布の異なる AuNCsを調製した。これらの試料に対し、1-Dodecanethiol 40 mM ヘキサン溶液を用いて、IL 中からの
AuNCs 抽出を試みた。今回用いた IL とヘキサンは室温において二相分離する。Au-IL 3.6ml に対しヘキサン溶液 1.8 ml を加え、激しく攪拌した後、24 h 毎にヘキサン相中の粒径を評価した。
粒径評価は小角 X 線散乱 (SAXS) 測定によって行い、あわせて UV-Vis 吸収測定も行った。
【結果と考察】攪拌後 24 h 程度で、ヘキサン相の着色が見られ、AuNCs のヘキサン相への移動を確認することができた。粒径の異なる 2 種類の試料について抽出を試みた結果を Fig. 1 に示す。ここから、粒径の違いにより、AuNCs の抽出速度に大きな差があることが分かる。このようにして得られた種々のサンプルについて SAXS 測定を行い、球状粒子を仮定した理論散乱曲線を用いてフィッ
ティングすることで粒径分布を導出した。抽出前の IL 中、抽出後のヘキサン相における粒径分布をFig. 1 抽出試料(41 倍希釈)Fig. 2 に示した。図を縦に貫く点線はイオン液体中におけるピーク位置を示している。どちらの試料においても、抽出時間 24 h での粒径分布は IL 中における粒径分布よりも小さいことが分かる。さらに抽出時間を増加させた場合の変化には、最大分布径2 nm と 4 nm の試料に違いが現れた。最大分布径が 2 nm 試料の場合には、48 時間までは元の粒径分布より小さな AuNCs が抽出されている。しかしながら、さらに長時間抽出して得られた試料では、はじめに IL 中で示したサイズよりも大きな粒子が得られている。一方最大分布径 4 nm の試料では抽出時間が増大しても粒径分布は大きく変化しない
ことがわかる。Fig. 1 より濃度が増大していることが明らかであるため、この試料の場合には 3 nm 程度の AuNCs が選択的に抽出されていることが考えられる。
これらの結果から、比較的簡単な操作によって、IL 相よりヘキサン相へと AuNCs の抽出が可能であることが明らかとなった。IL の静電的な保護能により、緩やかな保護を受けていたクラスターは、より強力な保護剤である 1-Dodecanethiol に保護されることで、ヘキサン相へ移動したものと考えられる。また、サイズが増大する場合もあるものの、粒子径によりヘキサン相への移動速度が異なり、これを利用した AuNCs のサイズ選別が可能であると考えられる。以上のように水－トルエン系と同様に、
IL－ヘキサン系においても、AuNCs のサイズ選択
的な抽出が可能であった。講演では、UV-Vis 測定の結果や、IL 中の AuNCs サイズを変化させた場合の測定結果も加えて、その抽出時間依存性について

Fig.2 粒径による抽出時間依存性議論する。
【参考文献】[1] T.Torimoto, K.Okazaki, T.Kiyama, K.Hirahara, N.Tanaka, S.Kuwabata, Appl. Phys. Lett.,
2006, 89, 243117.
[2] Y. Hatakeyama, M. Okamoto, T. Torimoto, S. Kuwabata, K. Nishikawa, J. Phys. Chem. C, 2009,

113, 3917.
[3] Y.Hatakeyama, S.Takahashi, K.Nishikawa, J. Phys. Chem. C, 2010, 114, 11098.
[4] Y. Hatakeyama, K. Onishi, K. Nishikawa, RSC Advances, 2011, 1, 1815.
[5] H. Yao, O. Momozawa, T. Hamatani, K. Kimura, Chem. Mater., 2001, 13, 4692

4P-078

共鳴フェムト秒誘導ラマン分光による光受容蛋白
ＰＹＰの励起状態構造ダイナミクスの研究(阪大ＶＢＬ1, 東北大院・理 2)

中村亮介 1, 濱田格雄 1, 阿部健太 2, 吉澤雅幸 2

Ultrafast structural evolution in the electronic excited state of
PYP studied by femtosecond stimulated Raman spectroscopy(Osaka Univ.1, Tohoku Univ.2) Ryosuke Nakamura1, Norio Hamada1, Kenta Abe2,
Masayuki Yoshizawa2【序】光受容蛋白質 Photoactive Yellow Protein (PYP)は、中A) WT

心色素である p-クマル酸（図１）が光励起されると、トラン

Glu46

O

Chromophore S

O H

ス-シス異性化を含む一連の構造変化がフェムト秒～サブ秒

O

O

H

の幅広い時間スケールで進行する。これまで、フェムト秒誘

N
H

O

Cys69

H

Tyr42

O

導ラマン散乱測定によって電子励起状態における共鳴ラマン

Thr50

信号を初めて観測し、光励起直後の構造変化について報告を行ってきた[1]。今回、p-クマル酸を保持している水素結合ネB) E46Q
Gln46

O

ットワークの機能的役割を明らかにすることを目的とし、特Hに p-クマル酸のフェノール部位と周囲のアミノ酸残基との水素結合に着目し、野生型（WT）および変異型(E46Q)との比

O

S
O

H

Tyr42

N
H

O

Cys69

H
O

Thr50

較を行った。
【実験】チタンサファイア再生増幅器（1kHz, 1mJ）からのChromophore

N H

図１ p-クマル酸と周囲のアミノ酸残基の構造。WT (A)と E46Q (B)。

パルス光を３分割して、励起光、ラマン励起光、プローブ光を発生させた。励起光（460 nm, 100fs）は、パラメトリック増幅（OPA）と和周波混合によって生成した。プローブ光は、サファイ
ア基板によって発生させた白色光を利用した。測定系の時間分解能は 150fs である。ラマン励起
光は狭帯域 OPA によって発生し、520 – 650 nm で波長可変とした。今回は主に中心波長 520.5 nm(幅 32 cm1)のラマン励起光を用いた。通常、第１励起状態 S1 の共鳴ラマン散乱を取得するには、
S1 – S2 遷移に共鳴したラマン励起光を用いる。われわれが用いたラマン励起光は誘導放出過程 S1– S0 遷移に共鳴している。この場合、ストークス側の信号は、逆ラマン過程によって励起されたS1 準位の振動を観測することになり、スペクトル形状は分散型になる。一方、アンチ・ストーク
ス側は、ラマン過程によって励起された S1 準位の振動を観測することになり、スペクトル形状は単純なローレンツ型になることが、これまでの研究から分かっている[1b]。従って、S1 準位の誘導ラマン散乱に関してはアンチ・ストークス側のスペクトルを議論する。

【結果と考察】図２Ａに、WT の基底状態 S0 における誘報告されている自発ラマン散乱スペクトルとほぼ一致している。1450 – 1600 cm1 には、芳香環およびエチレン部

WT
E46Q

(Stokes)Raman Gain

導ラマン散乱スペクトルを示す。この結果は、これまでにA Ground Stateの C-C, C=C 伸縮振動、1100 – 1350 cm1 には C=C（またB Excited State

0.1 ps(anti Stokes)

は C-C）と面内 CH rocking 振動が結合したモード、1000ていることが分かる。S0 準位における両者のスペクトルの

1555

1290

違いは、これまでの報告と一致している[3]。1555 cm1 の

1165

830

1.0 ps1055

として、E46Q の 1555 cm1 のラマン強度が大きく減少し0.3 ps990

1290 cm1 のモードで規格化して図示した。明らかな違い

900る[2]。同じ図に E46Q の誘導ラマン散乱スペクトルを

Raman Losscm1 以下には、skeleton モードや面外振動などが含まれ振動モードはチロシン Y19a に対応する p-クマル酸芳香環の伸縮モードであり、フェノール部位の水素結合状態に大きく依存することが知られている。
図２Ｂに励起状態 S1 における誘導ラマンスペクトルを

800

S1 では、WT と E46Q とがほぼ同じスペクトルとなってい1.6

I1555 /I1290

E46Q は S0 と S1 とであまり変化していない。結果として、

1.6

その減衰時間は時間分解能とほぼ同程度の 150 fs である。
一方で、E46Q は、励起後、相対強度比はほとんど変化せず 0.98 であり、S0 での値 1.15 に近い。
以上の結果から、WT では励起後、Glu46 と p-クマル酸
フェノール部位との水素結合が、時間分解能 150 fs よりも

I1555 /I1290を詳細に見てみると、図３に示すように両者で異なること減少し、1 ps 以上の時間領域では一定値 0.92 に達する。

A WT

1.2
1

で示すように 1.64 である。励起後、相対強度比は急激に

1600

1.4

る。しかし、1555 cm1 のラマン信号相対強度の時間変化が分かる。WT の場合、S0 における相対強度比は図中矢印

1200
1400
-1
Raman Shift (cm )図２ 基底状態(A)および励起状態(B)における誘導ラマン散乱スペクトル。

示す。WT は基底状態 S0 に比べて、1555 cm1 のラマン信号の相対強度が大きく減少していることが分かる。一方、

1000

B E46Q

1.4
1.2
1
0

0.5

1

1.5

Delay Time (ps)
図３ 1290 cm1 のラマン信号に対する 1555 cm1 の相対強度。WT (A)と E46Q (B)。図中の実線は 0.15psの時定数で減衰する指数関数。矢印は基底状態での相対強度比

速い時間で変化していると考えられる。これまでに報告されている中赤外フェムト秒過渡吸収分光[4]やピコ秒時間分解紫外共鳴ラマン散乱[5]との結果と合わせて議論する。
References[1] (a) Nakamura, R et al., Proceedings of The 15th International Conference on Time-Resolved
Vibrational Spectroscopy, 2011. (b) Nakamura, R et al., Proceedings of The 18th International
Conference on Ultrafast Phenomena, 2012.
[2] Unno, M.; Kumauchi, M.; Tokunaga, F.; Yamauchi, S. J. Phys. Chem. B 2007, 111, 2719-2726.
[3] Zhou, Y. et al., J. Phys. Chem. A 2001, 105, 5719-5726.
[4] Groot, M. L. et al., Biochemistry 2003, 42, 10054-10059.
[5] Mizuno, M.; Kamikubo, H.; Kataoka, M.; Mizutani Y. J. Phys.Chem. B 2011, 115, 9306-9310.

4P079分子動力学計算
分子動力学計算による動力学計算による芳香族炭化水素受容体による芳香族炭化水素受容体と芳香族炭化水素受容体とその共役因子
共役因子の
複合体の
安定構造の
その
共役因子
の複合体
の安定構造
の探索

(豊橋技術科学大学大学院 1 、理化学研究所計算科学機構 2、東芝研究開発センター3)○宮城慧 1、佐篠和哉 1、村田享士郎 1、伊藤聡 2、石原-菅野美津子 3、栗田典之 1,*

Molecular dynamics simulations for searching stable structure ofAhR complex with co-factor protein(Toyohashi University of Technology1, AICS RIKEN2, R&D Center Toshiba Corporation3)◯Satoshi Miyagi1, Kazuya Sashino1, Kyoshiro Murata1, Satoshi Itoh2, Mitsuko Ishihara-Sugano3 andNoriyuki Kurita1,*【はじめに】
はじめに】
これまでに、芳香族炭化水素受容体 (AhR: Aryl hydrocarbon receptor)は、外来異物が生体内に取り込まれた際に、それを特異的に結合し、高い親和性で認識し、その情報を核に伝え、代謝酵素の発現を誘導するタンパク質であることが明らかになっていた。さらに、近年の細胞実験[1-2]により、AhR が細胞分化にも重要な役割を担っていることが明らかになった。例えば、ナイーブ T
細胞が Th17 細胞と制御性 T 細胞へ分化をする際、AhR による転写調節が関わっていることが報告された[1]。その結果、自己免疫疾患やアレルギーといった細胞内外の環境によって調節される免疫系の制御機構に AhR が関与することが初めて示された。ナイーブ T 細胞の分化の方向性は、
2 種 類 の AhR リ ガ ン ド に よ っ て 変 化 す る [2] 。 そ の １ つ は 外 来 性 リ ガ ン ド TCDD
(2,3,7,8-tetrachlorodibenzo-p-dioxin) で あ り 、 も う 一 方 は 内 在 性 リ ガ ン ド FICZ (6-formylindolo[3,2-b]carbazole) である。現状では、AhR に結合するリガンドの違いにより、AhR がどのように T細胞の分化を分岐させるかは、原子レベルでは明らかになっていない。また、AhR とリガンド間の相互作用解析の際に必須である AhR の立体構造は未解明である。
我々は、これまでに、マウス AhR の LBD (Ligand binding domain)の構造、及び LBD に TCDD が結合した複合体の水中での安定構造を、分子モデリング手法と古典分子力学法を用いて求めた[3]。
また、分子モデリングで用いる鋳型構造を代え、AhR と様々なダイオキシン類を結合させた複合体の構造を水中で最適化し、実際の細胞実験と比較した結果も発表されている[4]。しかし、AhR
による異物認識機構においては、AhR は異物を結合した後、DNA の対象領域に結合するために共役因子と二量体を形成し、その二量体が DNA に特異的に結合すると考えられている。本研究では、
この機構の詳細を解明する目的で、AhR と共役因子である ARNT (Aryl hydrocarbon receptor nuclear
translocator)の二量体を分子モデリングによって作成し、リガンドの有無による二量体間の結合特性の変化を高精度分子シミュレーションにより解析した。
【計算手法
計算手法】
手法】
本研究では、リガンドとして外来性の TCDD と内在性の FICZ を採用し、ヒューマン AhR とARNT を含めた複合体の安定構造を探索し、AhR へのリガンド結合が、AhR と ARNT の二量体形
成にどのような影響を与えるかを解析した。AhR の LBD の立体構造は、タンパク質の構造予測プ
ログラム MODELLER を用い、ARNT に関しては、これまでに AhR の鋳型として用いていた構造の ARNT の部分構造をそのまま使用し、二量体の立体構造を作成した。更に、タンパク質リガン
ドドッキングプログラム Autodock を用い、二量体の様々な位置にリガンドをドッキングし、複合体の候補構造を多数作成し、古典分子力学計算プログラム AMBER9、及び Fragment molecularorbital (FMO)[5]法を用い、最安定な複合体の構造を決定した。その後、リガンド有無及び種類の変化により、二量体間の結合特性がどのように変化するかを明らかにするため、古典分子動力学(MD)計算を用い、複合体の構造変化をより広い領域で探索した。最後に、FMO 計算を用い、二量体形成に重要なアミノ酸の解明、リガンド有無による二量体間の結合特性の変化を解析した。
【結果と結果と考察】
考察】
Autodock を用いて求めた AhR-ARNT 二量体と TCDD の複合体、及び FICZ の複合体には、共に 4 個の代表構造があり、それぞれリガンドの結合位置が異なる。それらの構造を水中で最適化
し、FMO 計算によってエネルギーを計算し、Figures 1a, 1b に示す最安定構造を得た。従来のドッキングシミュレーション[4]により、AhR と TCDD の複合体においては、AhR の 7 つのアミノ酸が TCDD に相互作用することが明らかになっている。今回作成した Figure 1a の構造においても、同じ７つのアミノ酸が TCDD 周囲に存在し、作成した複合体構造の妥当性を確認できた。
これらの構造を初期構造にして、300K の水中で、MD 計算を 1nsec 行った結果を、Figures 2a, 2bに示す。Figure 1 と比べると、リガンド結合部位近傍の α ヘリックス構造が崩れている。また、二量体の AhR 及び ARNT の重心間距離は Figure 2a で 7.18 Å、Figure 2b で 4.88 Å であり、リガン
ドを含まない二量体（Figure 2c）では重心間距離は 13.75 Å であるので、リガンドの有無と種類により、
AhR と ARNT 間の相互作用が変化することが明らかになった。これらの構造に対し、FMO
計算により、AhR と ARNT 間の特異的相互作用を解析した結果は、当日のポスターで発表する。

(a) AhR-ARNT + TCDD(b) AhR-ARNT + FICZ

Figure 1 MM-optimized structures of the complexes with AhR-ARNT dimer and TCDD/FICZ(a) AhR-ARNT + TCDD(b) AhR-ARNT + FICZ(c) AhR-ARNT

Figure 2 MD-simulated structures of the complexes with AhR-ARNT dimer and TCDD/FICZ[1] F. J. Quintana, et al., Nature, 2008, 453, 65. [2] M. Veldhoen, et al., Nature, 2008, 453, 106.
[3] S. Miyagi, et al., Int. J. Quantum Chem., 2012, 112, 289.
[4] I. Motto, et al., J.Chem.Inf.Model, 2011, 51, 2868. [5] K. Kitaura, et al., Chem.Phys.Lett., 2001,336,163.

分子シミュレーション
分子シミュレーションによるシミュレーションによる芳香族炭化水素受容体による芳香族炭化水素受容体への芳香族炭化水素受容体へのリガンドの安定結合配置のリガンドの安定結合配置の探索

4P080(豊橋技術科学大学大学院 1、理化学研究所計算科学研究機構 2、東芝研究開発センター3 )村田享士郎 1、宮城慧 1、佐篠和哉 1、伊藤聡 2、石原-菅野美津子 3、栗田典之 1

Molecular simulations for searching stable binding-site ofvarious ligands into AhR protein(Toyohashi University of Technology 1, RIKEN AICS 2, TOSHIBA R&D Center 3)
Kyoshiro Murata1, Satoshi Miyagi1, Kazuya Sashino1, Satoshi Ito 2, Mitsuko Ishihara-Sugano3,
and Noriyuki Kurita 1【はじめに】
芳香族炭化水素受容体（AhR：Aryl hrydrocarbon receptor）は、リガンド依存性転写因子であり、外来異物や体内に存在する様々なリガンドを結合し、その情報を核に伝える。最近の研究により、
AhR は免疫細胞の分化や薬物代謝など様々な機構に関与している事が明らかになっている[1]。
AhR に 結 合 す る リ ガ ン ド に は 、 Fig.1 に 示 す 外 来 性 の ダ イ オ キ シ ン TCDD 、 及 び
FICZ(6-formylindolo[3,2-b]carbazole)、ITE(2-(1'H-indole-3'-carbonyl)-thiazole-4-carboxylic acid methyl
ester)などの生体内に存在する内因性のものがある。外来性リガンドは、毒性発現に関与する事が知られているが、内因性リガンドが関与する細胞の分化・再生機構及び免疫機構については未解
明な部分が多い。これらを解明するには、AhR の立体構造及び様々なリガンドとの特異的相互作
用を、原子・電子レベルで解明することが重要である。
我々は、これまでに、マウス AhR の LBD (Ligand binding domain)の構造、及び LBD に TCDD、
β-NF (β-Naphthoflavone)、FICZ、ITE が結合した 4 種類の複合体の水中での安定構造を、分子モデリ

ング手法と古典分子力場法を用いて求めた。さらに、LBD とリガンド間の相互作用の強さを、Ab
initio fragment molecular orbital (FMO)法[2]を用いて解析し、生化学実験で求めた各リガンドの AhR
転写反応への影響を、計算結果を基に電子レベルで説明しようとした[3]。本研究では、これまでに多くの実験研究が行われているラット AhR (以下 rAhR)へのリガンドの結合位置をより広く探索し、rAhR とリガンド間の特異的相互作用を明らかにした。
【計算手法】
本研究では、リガンドとして Figure 1 に示す 4 種類の化合物を採用し、rAhR の LBD とリガンド間の特異的相互作用を電子レベルで解析した。rAhR の LBD の立体構造は、アミノ酸並列情報を基に、タンパク質構造予測プログラム MODELLER を用い予測した。その際、鋳型構造として
HIF-2α (hypoxia inducible factors-2α)と人工リガンドが結合した PDB 構造(PDB ID: 3H82)を用い、人工リガンドを消去した構造に、タンパク質リガンドドッキングプログラム Autodock を用い、4 種類のリガンドをドッキングし、リガンドの安定結合位置を広範囲に探索した。さらに、幾つかの
rAhR+リガンド複合体の周囲に 8 Å に溶媒水を付加し、複合体の水和構造を古典分子動力学計算
プログラム AMBER9 を用い最適化した。最後に、FMO 法を用い、rAhR とリガンド間の特異的相互作用を、電子レベルで解析した。
【結果と考察】Autodock を用い、rAhR に TCDD をドッキングした結果を Fig. 2a に示す。生成した 256 個の候補構造は、TCDD がリガンドポケット内に入った構造(2 ケース)、ポケットの入り口付近に結合した構造、rAhR 外側の α-helix あるいは β-sheet 付近に結合した構造の 5 種類に分類できる。これらを水中で最適化した構造のエネルギーを FMO 計算で求め、最安定構造を決定した結果を Fig. 2bに示す。その構造の中で、TCDD がリガンドポケット内に入った構造が最安定になり、今回の計算手法の妥当性が検証できた。
さらに、得られた構造の妥当性を検証するため、Fig. 2b の構造を従来の結果と比較した。Motto
らの研究[4]では、rAhR と TCDD の結合に重要な 7 つのアミノ酸が明らかになっている。今回求めた rAhR+TCDD の最安定構造に対し、rAhR の各アミノ酸と TCDD 間の相互作用エネルギー(I.E.)を、FMO 法を用いて計算した結果、上記の 7 つのアミノ酸が TCDD と強い引力相互作用をしており、今回求めた構造の妥当性が確認できた。また、アミノ酸置換の実験[4]で、重要とされた他のアミノ酸に関しても、今回の計算で、TCDD と強く相互作用することが明らかになり、今回の一連の計算手法の妥当性が検証できた。
そこで、同じ計算手法を用い、rAhR に他の 3 種類のリガンドをドッキングし、rAhR との相互作用を調べた。まず、rAhR とリガンド間の結合の強さを見積るため、結合エネルギーを計算した。
その値は、30.8 (TCDD)、38.0 (β-NF)、51.0 (FICZ)、55.8 kcal/mol (ITE)であり、この結果は TCDDと FICZ 及び ITE の結合エネルギーの大小関係に関しては、従来の計算結果[5]と定性的に一致する。更に、rAhR とリガンド間の特異的相互作用を明らかにするため、rAhR の各アミノ酸とリガンド間の I.E.を FMO 計算により求めた。その結果、rAhR の Gln381 は、4 種類のリガンド全てと強く引力相互作用し、特に FICZ との相互作用が強い(-23 kcal/mol)ことが明らかになった。また、
Tyr320 は TCDD 及び ITE とのみ強く相互作用し、特に ITE と強く(-20 kcal/mol)相互作用する。従って、FICZ 及び ITE と rAhR 間の結合エネルギーが大きくなる主な原因は、Gln381 と FICZ、及び Tyr320 と ITE 間の特異的結合であることが、今回の計算で明らかになった。詳細は、当日のポスターで発表する。
【参考文献】[1] E. A. Stevens, et al., Nature, 2008, 453, 46.
[2] K. Kitaura, et al., Chem. Phys. Lett., 2001, 336, 163.
[3] K. Sashino, et al. , CBI/JSBi2011, 2011, Kobe, Japan.
[4] I. Motto, et al., J. Chem. Inf. Model, 2011, 51, 2868.
[5] H. Bisson, et al., J. Med. Chem., 2009, 52, 5635.

(a) TCDD(b) β-NF(c) FICZ(d) ITE
Figure 1 Structures of ligands for AhR.

(a)

(b)

Figure 2 (a) Various binding-sites of TCDD to LBD of rAhR produced by Autodock, (b) most stablestructure of solvated rAhR+TCDD optimized by AMBER9.

4P081BLUF ドメインの光活性化に伴う水素結合構造の赤外分光解析（名工大院工 1,名工大若手イノベータ 2,東邦大薬 3,光産業創成大 4）岩田達也 1,2,伊藤奨太 1,
伊関峰生 3,渡辺正勝 4,神取秀樹 1

FTIR study of hydrogen-bonding structures during the photoreactionof the BLUF domains(Grad. Sch. Eng., Nagoya Inst. Tech.1, Ctr. Fost. Yng. Innov. Nagoya Inst. Tech.2, Dept. Pharmacol.,Toho Univ.3, Grad. Sch. for Creation of Photonics Indust.4) Tatsuya Iwata1,2, Shota Ito1,
Mineo Iseki3, Masakatsu Watanabe4, Hideki Kandori1【序】BLUF（sensor of blue-light using FAD）ドメインは発色団として FAD（フラビンアデニ
ンジヌクレオチド）をもつ光受容体であり、バクテリアや真核生物において転写の抑制や光驚動反応など様々な機能に関わっている 1,2。BLUF ドメインは酸化型 FAD を暗状態（BLUF）に持ち、光励起されると電子･プロトン移動の後に再び酸化型 FAD をもった光反応中間体（BLUFred）になる 3。つまり、BLUF ドメインは他の光受容体とは異なり光活性化において発色団である FAD の化学構造が変化しないという特徴を持っている。BLUFred の生成に伴いFAD 近傍に位置する Tyr と Gln との水素結合構造が変化することが活性化に必須であると考えられているが、BLUF 及び BLUFred における水素結合ネットワーク構造は不明である（Figure 1）。そこで我々はフーリエ変換赤外（FTIR）差スペクトル分光法を用いた研究を開始した。特に注目したのは、水素結合ドナーとなる振動が現れる X-H 伸縮振動の領域であり、BLUFred において Tyr の O-H 伸縮振動が 2800-2600 cm-1 に現れることを発見した 4。
文献にないほどの低い振動数から、活性化によって Tyr が異常に強い水素結合を形成することを示している。本研究では、さらに研究を進め、BLUF 及び BLUFred における水素結合
構造の全貌を明らかにするため、系統的な同Figure 1 想 定 さ れ る BLUF ド メ イ ン のFAD, Tyr, Gln の構造位体標識を用いた赤外分光解析を行った。
【実験】紅色細菌由来の AppA-BLUF ドメインは His タグを N 末端側に融合させたタンパク質として大腸菌により発現させた。同位体標識試料の作製は 13C グルコース、15NH4Cl の入った培地を用いて行った。BLUF は封入体を生成したので、回収後これを変性させ、FAD 再生
buffer により再生後、Ni-NTA で精製した。分光測定においては、作製した乾燥フィルムに 1 μlの軽水または重水を添加し、水和フィルムの状態で実験を行った。試料を 260 K にセットし、青色光照射前後の赤外スペクトル変化を観測した。

【結果･考察】D2O 水和した試料に対して、BLUF（負の信号）と BLUFred（正の信号）における X-D 伸
縮振動の変化を Figure 2 に示す。Tyr-D4 標識試
料のスペクトルを非標識試料のものと比較したところ、BLUF で 2421 cm-1、BLUFred で 2029
cm-1 のバンドに同位体効果が見出されたことから、これらを Tyr 側鎖の O-D 伸縮振動と帰属した。対応する O-H 伸縮振動は BLUF で3292 cm-1 に現れることから 4、光活性化によって水素結合が異常に強くなる活性中心の Tyr
は、暗状態においても水素結合を形成していることがわかった。
次に Gln の C=O、C=N 伸縮振動の帰属を行った。以前、Gln が光反応に伴ってケト型（C=O）
からエノール型（C=N）に異性化が起こるのではないかという提案があった 5。そのため、アFigure 2

FTIR 差スペクトル X-D 伸縮振動領域

ポタンパク質のみを 13C 標識することで FAD
由来のシグナルとアポタンパク質由来のシグナルを分離した上でさらに、13C,15N 二重標識
試料を測定した（Figure 3）。その結果 1641
cm-1,1624 cm-1 のバンドがそれぞれ 5, 7 cm-1 低波数シフトした。このシフトは、Gln がエノール型に転移したことによる C=N 伸縮振動としては小さく、ケト型の C=O 伸縮振動が 15N 標
識の影響を間接的に受けたものであると解釈した。また、Gln の標識試料の測定からもエノール型の C=N 伸縮振動は観測できなかった。
よって光反応の前後で Gln はケト型の構造を

Figure 3 FTIR 差 ス ペ ク ト ル
C=O、C=N 伸縮振動領域

していると考えられる。さらにアポタンパク質だけを 15N 標識した試料と FAD も含めて 15N
標識した試料を比較することで FAD の ND 伸縮振動が BLUF で強い水素結合を形成しており、
BLUFred でさらに水素結合強度が強くなることがわかった。
本発表では 13C、15N、13C,15N、Tyr-D4、Gln-α-15N などの標識試料の赤外スペクトルをもとに活性中心の水素結合構造変化に関するモデルを提示したい。
【引用文献】
1. Han, Y., et al. Proc. Natl. Acad. Sci. USA 33, 12306-12311 (2004). 2. Yoshikawa, S., et al. Photochem.
Photobiol. Sci. 101, 727-731 (2005). 3. Gauden, M., et al.

Proc. Natl. Acad. Sci. USA 29, 10895-10900 (2006).

4. Iwata, T., et al J. Phys. Chem. Lett. 2, 1015-1019 (2011) . 5. Domratcheva, T., et al. Biophys. J. 94, 3872-3879(2008).

分子シミュレーション
βの分子シミュレーションによるシミュレーションによるアミロイド
によるアミロイドβ
アミロイドβタンパク質
タンパク質 AβC 末端変化
末端変化に相互作用変化のに伴う Aββ間の相互作用変化の解析

4P4P-082（豊橋技術科学大学大学院 1、村田製作所 2）
矢野篤志 1、岡本晃澄 1、野村和哉 1、檜貝信一 2、近藤孝志 2、神波誠治 2、栗田典之 1

Effects by addition of C-terminal residue on interactions among amyloid-β proteins:
Molecular simulation analyses(Toyohashi University of Technology1, Murata Manufacturing Co., Ltd.2)
Atsushi Yano1, Akisumi Okamoto1, Kazuya Nomura1, Shin’ichi Higai2, Takashi Kondo2,
Seiji Kamba2 and Noriyuki Kurita1【はじめに】
はじめに】
アルツハイマー病（AD）は、認知症の中で最も多い症状であり、その発症には、アミロイドβタンパク質（Aβ）の凝集が深く関係している。これまで、42 個のアミノ酸から成る、凝集性と毒性の高い Aβ42 が、その原因とされてきたが、近年、43 個のアミノ酸から成る Aβ43 について、
他の Aβより強い凝集性と毒性を示す事が報告されている [1]。本研究では、Aβ42 の C 末端にアミノ酸を付加した Aβ43 を作成し、古典分子力学（MM）及びフラグメント分子軌道（FMO）計算を用いて、Aβ間の相互作用の変化を解析、C 末端のアミノ酸付加が、Aβ間の結合にどのような影響を与えるかを明らかにした。
【計算手法】
計算手法】
本研究では、C 末端付加による Aβ間の相互作用変化を解析するため、Aβ42 五量体の実験構造（PDB ID: 2BEG）から Aβ42 二量体の構造を切り出し、その構造を基に、Aβ43 を含んだ二量体構造を作成した。Aβ42 五量体の実験構造から Aβ42 二量体を切り出す際、4 種類の隣接する Aβ42
二量体（A42B42、B42C42、C42D42、D42E42）を切り出し、各 Aβ42 二量体の周囲 8 Å に水分子を付加し、古典分子力場 AMBER99 を用い、水和構造を最適化した。最適化した Aβ42 二量体の各構造のエネルギーを FMO 計算（MP2/6-31G）により高精度に計算し、その結果に基づき、最安
定な Aβ42 二量体の構造（D42E42）を決定した。
次に、最安定な Aβ42 二量体構造 D42E42 に対し、D 鎖、E 鎖の C 末端に Thr を付加し、Aβ43を含んだ二量体の初期構造を作成した。その際、Thr を D 鎖のみに付加した構造を D43E42、E 鎖のみに付加した構造を D42E43、両方に付加した構造を D43E43 とし、Thr の側鎖の向きを 2 方向
考慮し、より安定な向きを決定した。これらの構造の周囲 8 Å に水分子を付加し、AMBER99 を用いて水和構造を最適化し、各構造の電子状態を FMO 計算により高精度に解析、Aβ二量体間の特異的相互作用が、Aβ43 の C 末端 Thr によりどのように変化するかを、電子レベルで解析した。
【計算結果
計算結果と
結果と考察】
考察】
水中で最適化した 4 種類の Aβ42 二量体の中で、最安定な構造（D42E42）を Figure 1 (a) に示す。
この構造は、他の 3 種類の構造と比較して、261.1 kcal/mol 安定である。この安定化の原因は以下のように考えられる。初期構造作成に用いた Aβ42 五量体の実験構造において、A 鎖、E 鎖は表面に存在し、周囲の溶媒との間で安定な構造を形成している。また、E 鎖中の Asp23 と Lys28 の側鎖間に水素結合が形成され、E 鎖は安定化しているが、A 鎖ではこの水素結合が存在しない。そのため、安定な E 鎖とそれに隣接した D 鎖を用いた二量体 D42E42 が、溶媒中では最安定になったと考えられる。現在、古典分子動力学（MD）シミュレーションを実行し、二量体の安定構造を

より広範囲に探索しようとしている。
この D42E42 最適化構造の C 末端に Thr を付加し、水中で最適化した構造を、
Figures 1 (b)−(d) に示す。Thr を付加した各構造と D42E42 との共通部分の RMSD を比較した結果、二量体全体では2.8 Å の構造変化が生じ、特に Thr43 を付加した D 鎖の C 末端付近で最大約 6.6 Å の構造変化が確認された。この原因は、一方の鎖に付加した Thr43 が、別の鎖の C 末端付近のアミノ酸側鎖と相互作用したためであると考えられる。
次に、Aβ間の結合の強さを予測するため、Figure 1 に示す最適化構造の電子状態を FMO 計算により計算し、Aβ間の結合エネルギーを求めた結果を Table 1 に示す。結合エネルギーは、D42E42
二量体で最小であり、D42E43、D43E42、D43E43 の順で大きくなる。D42E43 と D43E42 の結合エネルギーが異なる原因は、付加した Thr43 の側鎖の向きが関係している。より安定な D43E42 中
では、Figure 1 (b) に示すように、D 鎖に付加した Thr43 の側鎖は E 鎖の方向を向き、E 鎖のアミノ酸側鎖と相互作用するが、D42E43 ではそのような相互作用は存在しない。また、D 及び E 鎖に
Thr43 を付加した D43E43 では、両方の Thr43 が別の鎖のアミノ酸側鎖と相互作用し、より強い結合が実現している。さらに、Aβ中のどのアミノ酸が二量体形成に重要であるかを明らかにするため、D 及び E 鎖の各アミノ酸間の相互作用を解析した結果、D42E42 においては、D 鎖、E 鎖それ
ぞれの Asp23 と Lys28 が強く引力相互作用し、一方、D43E43 においては、Aβ中の Asp23 と Lys28
間の引力相互作用に加え、Thr43 同士が強く引力相互作用することが分かった。各アミノ酸間の相互作用の詳細、及び MD 計算により得た結果に関しては、当日のポスターで発表する。
Side chain of Thr43
D chain (D42)

D chain (D43)

E chain (E42)
E chain (E42)(a) D42E42
D chain (D42)

E chain (E43)(c) D42E43(b) D43E42

D chain (D43)

E chain (E43)(d) D43E43

Figure 1 Solvated structures of Aβ dimers optimized by AMBER-MM method
Table 1 Total energies of Aβ dimers and binding energies between Aβ

Binding energy
Total energy (kcal/mol)Complex D chain + water E chain + water
Water(kcal/mol)
D42E42 −44319378.4 −38547020.7 −38547091.4
261.9−32774995.6
D42E43 −44545379.4 −38547149.1 −38773237.8
293.7−32775301.3
D43E42 −44545319.4 −38773197.7 −38547078.1
300.5−32775256.8−32775233.4
310.1
D43E43 −44771492.0 −38773163.9 −38773251.4

Structure【参考文献】
参考文献】[1] T. Saito, et al, Nature, 14 (2011) 1023.

4P-083
第一原理分子軌道計算による第一原理分子軌道計算による
DNA 塩基対の異性化反応の遷移状態の塩基対の異性化反応の遷移状態の探索（豊橋技術科学大学大学院 1、ウクライナ科学アカデミー2）
野村和哉 1、干場靖浩 1、Victor I. Danilov2、栗田典之 1

Ab initio MO study on transition states
between wobble and tautomeric DNA base pairs(Toyohashi University of Technology1, National Academy of Sciences of Ukraine2)
Kazuya Nomura1, Yasuhiro Hoshiba1, Victor I. Danilov2 and Noriyuki Kurita1【はじめに】
はじめに】生命の遺伝情報を担う DNA 二重鎖の形成において、Watson と Crick が提案した標準的な塩基対（G-C, A-T）が、相補的な水素結合相互作用により正確に形成されることが重要である[1]。DNA
の複製過程において、新しい DNA 鎖が正しく合成されるためには、この塩基の相補性が必須である。しかし、DNA ポリメラーゼが、間違ったヌクレオチドを挿入することで、突然変異を引き起こすことがある。この原因の一つとして塩基対の互変異性による影響が考えられる。実際、DNA複製の際に、互変異性によって引き起こされるエラーの頻度は 10-4～10-5 であり、互変異性を引き起こす機構を解明することは非常に重要である。また、特定の塩基対がゆらぎ位置で出現する原因を理解するため、様々なコドン・アンチコドン相互作用の組み合わせが調査され[2]、12 量体のd(CGCGAATTTGCG)2 二重鎖において、ゆらぎ塩基対(wG-T)が 2 つ存在することが実証された[3]。
wG-T 塩基対は 2 つのプロトン移動により、enol 型 G（G*）あるいは enol 型 T（T*）を含む G-T
塩基対へ互変異性化することができる。
本研究では、wG-T 塩基対の互変異性の可能性を解明するため、第一原理分子軌道計算を用い、
wG-T から G*-T あるいは G-T*の塩基対に互変異性化する反応機構を解析した。
また、変異原性の強い塩基類似体の一つである Bromouracil（BrU）は、enol 型の互変異性体となり、A ではなく G と誤った塩基対を形成することが知られている。そこで、BrU の影響を明らかにするため、wG-BrU 塩基対に対し同様の解析を行い、BrU の異性化反応への影響を調べた。
【計算手法
計算手法】
手法】
本研究では、まず、wG-T、wG-BrU とそれらの互変異性化構造 G-T*、G*-T、G-BrU*及び G*-BrUの 構 造 を 第 一 原 理 分 子 軌 道 計 算 プ ロ グ ラ ム Gaussian09(G09) の M06/6-311++G(d,P) 法 、 及 び
MP2/6-31G(d,p)法を用い、真空中で最適化した。次に、G09 の QST2、QST3 法を用い、これらの塩基対構造間の遷移状態(TS)を探索した。QST3 計算では、反応前後の構造に TS の候補構造を付け加えた 3 つの構造を基に、TS 構造を探索した。その後、得られた TS 構造が遷移状態であることを確かめるため、振動解析を行い、虚数振動数が 1 つのみであることを確認した。最後に、TS構造を初期構造として、G09 の IRC 計算（ステップ数 30）を行い、TS 構造から反応前後の構造への構造変化を解析し、得られた TS 構造が 2 つの構造を結ぶ遷移状態であることを確認した。
【計算結果と計算結果と考察】考察】G09 の MP2/6-31G(d,p)法と QST 法により求めた 3 種類の G-T 塩基対間の遷移状態を Figure 1 に示す。求めた構造は、1 つの虚数振動数を持ち遷移構造である。M06/6-311++G(d,P)法を用いた場

合も、ほぼ同じ結果が得られている。TS(wG-T→G*-T)及び、TS(wG-T→G-T*)については、T のプロトンが G に転移し、T の構造が変化した後、G からのプロトン転移によって G*-T、G-T*の構造へと遷移するため、G にプロトンが 3 つ存在する構造が遷移状態となる。また、TS(G-T*→G*-T)については、プロトンの移動によって G-T*から G*-T の構造に変化するため、2 つの構造間のほぼ等距離にプロトンが存在している構造が遷移状態となる。
各構造の相対的なエネルギーを、Figure 2 に示す。黒字の値は MP2/6-31G(d,p)法で求めた結果、青字の値は M06/6-311++G(d,P)法で求めた結果である。MP2/6-31G(d,p)法による結果では、3 つの塩基対の中で、G*-T が最安定であり、wG-T、G-T*はそれぞれ、1.0, 2.3 kcal/mol 不安定である。
wG-T から G*-T 及び G-T*に変化する反応障壁が等しく 19.4 kcal/mol であり、G-T*と G*-T 間の反応障壁は 4.8 kcal/mol である。この結果から、wG-T から G*-T 及び G-T*への反応は同じ遷移構造
を通して起こり、反応確率も等しくなるが、最終的には、G*-T の構造がより安定であるため、G-T*
から G*-T への互変異性化が起り、G*-T の構造の割合が多くなる。M06/6-311++G(d,P)法で求めた結果も、同様の傾向にある。
また、wG-BrU、G-BrU*、G*-BrU の塩基対に関しては、wG-BrU から G-BrU*及び G*-BrU への反応障壁は 17.0 kcal/mol となり、wG-T 塩基対の場合と比べ 2.4 kcal/mol 低くなる。従って、Br の影響により反応障壁が下がり、wG-T よりも異性化し易くし、enol 型塩基を生成し易くすることが明らかになった。その結果、BrU の影響で突然変異が起こり易くなると考えられる。計算結果の詳細は、当日のポスターにて発表する。

1.47Å

2.10 Å

2.10 Å

1.90

1.90
2.02

2.02

1.94

1.94

1.37

1.81

1.86

1.86

(a) TS(wG-T→G*-T)(b) TS(wG-T→G-T*)(c) TS(G-T*→G*-T)

Figure 1 Structures of transition states obtained by QST and MP2/6-31G(d,p) methods in G09
TS (wG-T → G-T*)

TS (wG-T → G*-T)

TS(G-T*→→G*-T)

G*-T

19.3
16.6

wG-T

19.3
16.6

G-T*

G-T*

G*-T
4.8
5.3

-2.3
-1.6

0.0 kcal/mol

0.1
-0.6

0.1
-0.6

-2.3
-1.6

Figure 2 Schematic energy diagrams for the transitions between wG-T, G-T* and G*-T: relative energies
obtained by MP2/6-31G(d,p) and M06/6-311++G(d,P) are shown in black and blue, respectively.
【参考文献】
参考文献】[1] F. H. C. Crick, Proc. R. Soc. B167 (1969) 331-347.
[2] M. D. Topal, J. R. Fresco, Nature 263 (1976) 289-293.
[3] T. A. Early, J. Olmsted III, D. R. Keams, A. G. Lezius, Nucleic Acids Research 5 (1978) 1955-1970.

4P-084
変性したプロテイン A の構造解析：CW とパルス ESR 法（東北大多元研 1)、東大院 2)）阿部 淳 1)、新井 宗仁 2)、高橋 聡 1)、大庭 裕範 1)、山内 清語 1)

A Structural Analysis of Denatured Protein A: CW- and Pulsed ESR（IMRAM, Tohoku Univ. 1)、Tokyo Univ. 2)）Jun. Abe 1)、Munehito Arai 2)、Satoshi Takahashi 1)、Yasunori
Ohba 1)、Seigo Yamauchi 1)【序】タンパク質フォールディングの分子機構を理解することで、未知タンパク質の構造予測やフォールディング病の対策など多くの応用が期待できる。アミノ酸 100 残基以下の単一ドメインからなるタンパク質のフォールディングは、天然状態と変性状態の二状態の平衡で示される（図 1）。ここで、変性状態におけるタンパク質は、ランダムコイルの性質を持つと考えられていた[1]。しかし、

天然状態

A55C

変性状態

A55C

N22C

近年になって変性状態における残余構造の存在を示唆する報告[2, 3]が多くなされるなど、ランダムコイル仮説に対する疑N22C

問も投げかけられている。このように変性状態の構造の解明はフォールディング

図1

プロテイン A の天然状態と変性状態の二状態平衡の出発点を理解し、変性したタンパク質の精密なモデル化のために非常に重要である。本研究ではフォールディング機構が良く知られているプロテイン A の B ドメイン(BDPA) [4-6]を対象として、CW-ESR とパルス ESR（二量子コヒーレンス ESR）を用いて変性状態の構造研究を行った。
【実験】 BDPA の 1 ヶ所にシステイン変異を導入した Y15F/A55C 置換体と 2 か所にシステ
イン変異を導入した Y15F/N22C/A55C 置換体にスピンラベル剤を反応させて、シングルラベ
ル体（A55C）とダブルラベル体（N22C/A55C）を作成した（図 1 に示した部位にラベル化)。
このときシングルラベル体の BDPA 濃度は 15~35 μM、ダブルラベル体の濃度は 42~70 μM と
Bradford 法より定量した。変性剤濃度 0 M のタンパク質 6は 150 mM の NaCl を含む pH 7.0 の MOPS 緩衝溶液で測 5 5 M GdmCl
定し、変性タンパク質の測定には 0.5~5 M の変性剤（塩酸
グアニジウム塩、GdmCl）を加えた。変性剤を含んだ緩衝溶液への交換は Sephadex G-25 desalting column を用い

4 M GdmCl
4

3 M GdmCl
3てゲルろ過により交換した。CW-ESR スペクトルの測定 2

2 M GdmClに は 、 JEOL の X バ ン ド (9.5 GHz) CW-ESR 装 置 1

1 M GdmCl（JES-FE2XGS）を用いて室温で測定した。また二量子コ

0 M GdmCl0

ヒーレンス ESR の測定は、JEOL と共同で開発した Ku バ
ンド(17.5 GHz) パルス ESR 装置を使用し、6 パルス系列

3240

で 64 ステップの位相サイクルを用いて 50 K で行った。

図 2【結果と考察】 室温でのシングルラベル体の CW-ESR 測

変性剤濃度依存性

3260

3280

3300

3320

Magnetic Field (G)CW-ESR スペクトルの

定結果を図 2 に示す。CW-ESR の線幅解析から、GdmCl 0~2 M に比べて 4~5 M の線幅が小さ
くなり、この原因はタンパク質主鎖の運動性が大きくなったことに起因していると考えられる。また GdmCl 4~5 M の範囲で線幅に変化がなかった。スペクトルの線幅解析より変性剤濃1.0

光の結果[5]と良い一致を示した。

0.8

ダブルラベル体での二量子コヒーレンス ESR（測定
温度 50 K）の結果から GdmCl 0 M では、わずかに双極Echo Amplitude度に対して線幅をプロットした結果は既報の円二色偏a0M GdmCl

0.6

0.4
0.2

子相互作用による振動が見られた（図 3a）。一方で0.0

GdmCl 1~5 M では（図 3c）明瞭な振動が見られなかっ1.0

た。これは、幅広い距離分布を持つことを示している。

0.8

また GdmCl 1~5 M の間でほとんど変化がみられなかった。GdmCl 0.5 M では（図 3b）早い成分と緩やかに減衰する成分の重ね合わせになっている。このことから
GdmCl 0.5 M では GdmCl 0 M と 1~5 M での双極子変調

Echo Amplitude

0.00

0.05

0.10

0.15

0.20

Time (µs)

b0.5M GdmCl

0.6

0.4
0.2

0.0
0.0

0.2

0.4

0.6

0.8

1.0

Time (µs)

cの重ね合わせになっていると示唆される。二量子コヒー
レンス ESR での変性中点は 0~1 M の間であり、既報の円二色偏光より得られた変性中点[4]に比べて小さくなっている。
この原因は二量子コヒーレンス ESR では 50 K
で測定を行うため試料を凍結する必要があり、それにより低温変性した可能性、あるいは凍結することにより変性 剤 濃度 が局 所的 に増加 し た可 能性 が推 測され る 。
GdmCl 0 M で距離分布解析を行うと二つの距離成分が現われ、2 つのガウス関数のフィッティングからそれぞれ平均距離が 1.51、2.18 nm であり、距離分布は 1.00、0.41 nm

図3

二量子コヒーレンス ESR の測定結果；a. 変性剤濃度 0 M； b.
変性剤濃度 0.5 M；c. 変性剤濃度1~5 Mであった。既報の NMR 構造[6]で主鎖の構造が異なる構造を基にしてスピンラベル側鎖だけを動かして計算した値との比較から、タンパク質の主鎖の構造の違いにより距離分布が二成分になると示された。GdmCl 5 M の実測と比較するために、ランダムコイルを仮定した距離分布での双極子変調を計算した。このとき計算した双極子変調は擬似永年項を入れた一次元の二量子コヒーレンス ESR のエコー信号である[7]。ランダムコイルモデルは理想鎖と非理想鎖の二つを考えた。理想鎖、非理想鎖を仮定した場合、共に実測よりも早い減衰となり、実測よりも短い距離の成分が多いことが示された。このことから GdmCl 5 M では理想鎖や非理想鎖とは異なる広がった構造を取っていることが示唆される。また GdmCl 0 M の距離分布関数が重なる領域が小さいため、変性状態は天然状態とは異なる構造であると考えられる。
GdmCl 5 M と GdmCl 1~4 M では二量子コヒーレンス ESR の結果が同じであることから、この濃度領域で同様の構造を取っていることが推測される。
【参考文献】[1] C. Tanford et al (1966) JBC, 241, 1921, [2] D. Shortle et al (2001) Science, 293,487, [3] E.
Sherman and G. Haran (2006) PNAS, 103, 11539, [4] F. Huang et al (2009) Biochem. 48, 3468, [5] S. Sato et al(2006) JMB, 360, 850, [6] S. Sato et al (2004) PNAS, 101, 6952, [7] S. K. Misra et al (2009) Appl. Magn. Reson.
36, 237,

4P085Elongation 局所構造最適化による全系の効率的安定構造探索(九大院・総理工1, SCNU2, JST-CREST3) Liu Kai1, 折本 裕一1, Gu Feng Long2, 3, 青木 百合子1, 3

Efficient search for optimized geometry of whole system by elongation local
optimization method(Kyushu Univ.1, SCNU2, JST-CREST3) Kai Liu1，Yuuichi Orimoto1，Feng Long Gu2，Yuriko Aoki1,3

I. Introduction
One of the most successful applications of molecular quantum mechanics is the reproduction and prediction ofmolecular conformation. However, it’s still a grand challenge to perform the geometry optimization on large system
in practices. Due to a large number of atoms, the standard ab initio calculations are beyond the computational reach.
Currently, some fragment-based methods have been successfully proposed and developed to calculate large system
within a reasonable cost, for example: divide and conquer (DC)1, elongation method (ELG)2, fragment molecularorbitals (FMO)3 and systematic fragmentation method (SFM)4. Despite the complicate interactions of large system,
it’s very difficult to locate the equilibrium structure. Here we present the implementation of geometry optimization
based on the elongation method, which is famous for its high accuracy and efficiency.

II. Theoretical approach
A. Elongation method
ELG procedure generally is analogous to experimental polymer
chain synthesis as shown in Fig.1. First, a suitable size of initialmonomers (starting cluster) is chosen to initiate the ELG procedure.
The canonical molecular orbitals (CMOs) of the starting cluster,generated by conventional SCF calculation, are transformed to an
orthogonal atomic basis (OAO). In the following step, the OAO-based
density matrix (DOAO) is to partition into frozen region (A) and active
region (B). The frozen region is assumed to be far away from thechain propagation point, while the active region is consisted of theremaining part of the starting cluster. After the separate
diagonalization of the subspace DOAO (A) and DOAO (B), a set of
regional localized molecular orbitals (RLMOs) for both A and B
regions are obtained. Then, one attacking unit (C) is added to theFig.1 The flowchart of elongation method
chain propagation point for the next ELG step. Because of thenegligible interactions between A and C, only B and C will be included in the ELG SCF calculations. After the ELG
HF-SCF converges, the CMOs of B and C regions will be localized again to form a new frozen region (RLMOs B’)and a new active region (RLMOs C’). Then a new attacking unit (D) is added to repeat the above procedures until
the desired length is reached. The important feature of the ELG method is that the equations during SCF are solved
only for small subunits instead of the whole system, and the calculations of two-electron repulsion integrals (ERIs)
between A and M can be partly omitted by cutoff technique.

B. Elongation geometry optimization
In the framework of HF calculation, the first derivative (gradient) of total energy (E) with respect to the nuclear
coordinate X A at atomic orbitals basis, can be written as:𝝏𝑯𝒄𝒐𝒓𝒆𝝏𝑺𝝁𝝂𝝏(𝝁𝝂||𝝈𝝀) 𝝏𝑽𝑵𝑵
𝝏𝑬
𝟏𝝁𝝂= � 𝑫𝝂𝝁+ � 𝑫𝝂𝝁 𝑫𝝀𝝈+− � 𝑸𝝂𝝁
𝝏𝑿𝑨
𝝏𝑿𝑨
𝝏𝑿𝑨
𝝏𝑿𝑨
𝝏𝑿𝑨
𝟐
𝝁𝝂

𝝁𝝂𝝀𝝈

𝝁𝝂

where (µν||σλ) denotes two-electron integral, D is the density matrix, Hcore and S correspond to the core Hamiltonianand overlap matrices of system, respectively. The nuclear-nuclear repulsion is defined by V NN . The energy-weighted

N/2
density matrix Q of conventional method is defined asQ νμ = ∑i ni εi Cμi Cνi , where n is the matrix of occupancy
number (density matrix in MO representation). In the ELG method, the equations during SCF are solved only for
active subspaces defined by B and M fragments. Therefore, the non-diagonal εLMO
is employed in the ELG -OPT
ij
AO
LMO
method. As the coefficients of the whole system CABM consist of CA (coefficients of frozen part, localized by the
AOMO
ELG localization procedure) and CBM(coefficients of active region and attacking unit, transformed from CBM
after
the ELG SCF calculation). Therefore, the energy-weighted density matrix Q of the ELG-HF-OPT method in atomic
N/2
basis can be rewritten asQ νμ = ∑ij ni εij Cμi Cνj . Because of the boundary effects between A and BM arisen by thetails after the ELG localization procedure, the coordinates of the one unit of the BM region, which is the closest to
the frozen region, is fixed in the gradient calculation to reduce these effects.

III. Results and discussion
A. Non-bonding model system: (HF) n=48
The linear poly-hydrogen fluoride (poly-HF) molecules are optimized by both elongation and conventional
restricted Hartree-Fork method (RHF) with different basis sets. The energy differences ∆E (∆E=Eelongtion-Econventional)of STO-3G, 6-31G and 6-31G(d,p) basis sets are -4.49×10-7, 9.09×10-9 and -3.26×10-5 Hartree/atom, respectively.
The negative values means the ELG-OPT locates an even lower ground state than conventional results. It indicatesthat for a flat energy potential surface of system, like linear poly-HF molecules, ELG-OPT may produce a more
promising candidate for the most stable geometry.

F

RF...H

H

F

RF...HH

6.00E-02
4.00E-02
2.00E-02

STO-3G
6-31G
6-31G(d,p)

CONV provides ametastable structure.

0.00E+00

0-2.00E-02

0-4,800.2

Total energy

8.00E-02

Final step of optimization iteration(au)-4,800.0

RF...H bond length difference(Å)
1.00E-01

20

40

60

6-31G(d,p)-4,800.4-4,800.680100

ELG-OPT
CONV-OPT

Eelg lower than Econv
Telg shorter than Tconv-4,800.8-4,801.0

5 10 15 20 25 30 35 40 45 50-4,801.2

E: Energy T: Time

Fig. 2 The differences of bond length and optimization iteration of poly-HF systems. ∆R=R

elongation

B. Bonding system: extended ployalanine.
As a model of large biosystem, 20 units ofextended alanines, is optimized by both the ELGHF-OPT and conventional method using 6-31G
basis set. The RMSD between the ELG andconventional optimized structure (excludinghydrogen) is 0.18. The bond length is also
compared in details and shown in Fig.3. The
difference of total energy between ELG-OPT andconventional result is 2.37×10-7 Hartree/atom.
These small differences in structure and totalenergy show the well reproduction of elongationgeometry optimization method.

-Rconventional

The difference of bond length(Å)
8.00E-04
6.00E-04

6-31G

4.00E-04
2.00E-04
5.00E-19-2.00E-04 0

20

40

60

80

100

-4.00E-04

Fig. 3 The difference of bond length between ELG andconventional optimized structure.

Compared to conventional results, it indicates that
the ELG-OPT can well reproduce the calculations and may locate a more stable structure than conventional one.
1. W. T. Yang, Phys. Rev. Lett. 1991, 66, 1438.
2. A. Imamura, Y. Aoki, K. Maekawa, J. Chem. Phys. 1991, 95, 5419.
3. K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayasi, Chem. Phys. Lett. 1999, 313, 701.
4. V. Deev, M. A. Collins, J. Chem. Phys. 2005, 122, 154102.

4P086

Elongation dynamics と反応部位局所遷移状態探索(九大院・総理工 1，SCNU2，JST-CREST3）Xie Peng1，折本 裕一 1，Liu Kai1，Yan Yun-an1，
Gu Feng Long2, 3，青木 百合子 1, 3

Elongation dynamics and local transition state search on
reaction site
(KYUSHU Univ.1, SCNU2, JST-CREST3) Peng Xie1，Yuuichi Orimoto1，Kai Liu1，Yun-an
Yan1，Feng Long Gu2, 3，Yuriko Aoki1, 3[Introduction]
Many essential biological processes are often accompanied by electron-related events suchas chemical reactions. The mechanisms of these chemical reactions are significant for theunderstanding of biomolecular functions. Ab initio molecular dynamics (AIMD) method is a
powerful tool to describe chemical reactions involving bond breaking and forming in whichelectrons play an important role, but the number of atoms handled with quantum
mechanics is much smaller than the number of atoms in a biological system for thelimitation of computer technologies. Elongation method (ELG)1,2 can efficiently calculate theelectronic structure of large aperiodic polymers with high accuracy on ab-initio level, so we
proposed a new method named Elongation molecular dynamics (ELG-MD) by combine theelongation method with dynamics method.
This ELG-MD method makes it possible toefficiently analyze the mechanisms ofchemical reactions for biomolecule. To
describe a chemical reaction, we may propose
several kinds of chemical mechanism. It is
important to make clear the transition statesalong the reaction pathway. To find andconfirm the saddle point, we need calculate
the Hessian matrix based on quantummechanical principles. For this purpose, wedeveloped ELG-Hessian algorithm which canefficiently calculate the local Hessian matrix
for the reaction site.
[Methodology]
In ELG-MD method, we use gear predictor
corrector (GPC) method to solve the equations Figure 1: The workflow of ELG-MDof motion for wavefunction and coordinate
dynamic variables within the ELG method. The procedure is as follows (Fig. 1): The whole
molecule will be divided into a number of units. The first step, we select the 1st unit (Aregion) and a suitable number of units (B region, to ensure A region has no interaction with
M region) as the starting cluster to perform AIMD calculation for a certain number of steps.

The second step, A region will be frozen and a new unit (M region) will be attached at the
tail site of B region. The local force on each atom of the active region (B+M region) will beobtained by ELG method, then the acceleration and new position of each atom will becalculated by GPC method. After a certain number steps of dynamics, we can repeat thesecond step till we obtain the whole system in equilibrium.
Hessian matrix calculation is considerable expensive for biological system when aquantum mechanical description is used, but the chemical reaction region is always on asmall part, the most regions do not interact with the reaction region. So in ELG-Hessian
calculation, we take the chemical reaction region as the last unit (N) of elongation. From thestarting cluster to the second last step, we only need to obtain local molecular orbitas for thefrozen region (A region). For the last step of ELG-Hessian,
we use full molecular orbitals to calculate the Hessian
matrix element for B+M region.

Hx y =∂2 Etotal∂x∂yx ∈ B + M,y ∈ B+M(1)[Results and Discussion]
Polyglycine with water was selected as a target system
to test the ELG-MD method. In each ELG step, we
performed 4000 steps dynamics calculation at HF/STO-3G
level. The initial structure and the final structure areshown in Fig. 2 (a) and (b) respectively. Fig. 2(c) showed
the structure obtained by ELG-optimization at HF/STOFig. 2: (a) The initial structure of poly3G level. By comparing the Fig. 2 (b) and (c), we foundglycine with water. (b) The structure
obtained by ELG-MD. (c) The structure
that: the structure in Fig. 2 (b) has more helix caused byobtained by ELG-optimization.
hydrogen bond forming than Fig. 2 (c). That means theprevious structure is more close to the real structure of polyglycine.
Before performing ELG-MD calculation, we need use ELG-Hessian method to find thetransition state of the reaction site. To examine the accuracy and efficiency of ELG-Hessian,
we compared the vibration modes of linear polyglycine (20 units) at HF/STO-3G level withconventional (CNV) Hessian calculation.
Table 1 The last five intra-molecular vibration modes for
The comparisons are shown in Table 1.
linear polyglycine Frequency unit (wavenumber)
ELG 4054.98 4055.27 4055.62 4056.32 4057.51 We can see the error of frequency is
CNV 4058.98 4060.24 4058.21 4058.85 4060.86 considerable small. The mostly timeconsuming part of Hessian calculation is CPHF part, we checked the CPU time for CPHF
calculation. We found that the CPU time for CPHF calculation can significantly drop, from615.5 second (CNV) to 389.8 second (ELG).
[Future plan]
ELG-Hessian and ELG-MD can work well now, as the next step, elongation transition
states (ELG-TS) method will be programmed for efficient investigation on chemical reaction.
1. Imamura, A.; Aoki, Y.; Maekawa, K. Journal of Chemical Physics 1991, 95, 5419-5431.
2. Gu, F. L.; Aoki, Y.; Korchowiec, J.; Imamura, A.; Kirtman, B. Journal of Chemical
Physics 2004, 121, 10385-10391.

4P087非局所励起状態に対する分割統治(DC)型線形スケーリング法の数値検証（早大先進理工1，早大理工研2，JST-CREST3）
野中 佑太郎1，吉川 武司1，中井 浩巳1,2,3

Numerical assesments of divide-and-conquer based linear-scaling method
for non-local excited states(Advanced Science and Engineering, Waseda Univ.1, RISE, Waseda Univ.2, CREST, JST Agency3)
Yutaro Nonaka1, Takeshi Yoshikawa1, Hiromi Nakai1,2,3【緒言】
現在までにフラグメント分子軌道法[1]や分割統治(DC)法[2,3]等の様々な分割型計算理論が開発されてきた。これらの手法を用いることにより、数千原子を超える大規模分子に対して高精度な量子化学計算を行うことが可能となりつつある。特に、基底状態においては分光学的精度を保証す
る CCSD(T)にまで対応させることに成功した[4,5]。しかし、励起状態に対する計算は未だ発展途上であり、現状では Charge Transfer (CT)型励起のような非局所励起の取り扱いは非常に困難である。
さらに、DC 法においては結合の切断をあらわに取り入れているため非物理的な励起が現れるという問題点も存在する。本発表では、励起状態計算の現状を打開する手立てとして、周波数依存分極率を求める手法である DC-TDCPHF/DFT[6]を利用して、その極から励起エネルギーを求めることができることを数値的に検証する。

【DC-TDCPDFT の理論】
DC法は全系をいくつかの部分系に分けて計算を行うことで計算コストを削減する方法である。
部分系の周りのバッファ領域を含めた局在化領域で部分系sの軌道を構築することで、周囲の効果を取り込むことができる。
DC-TDCPDFT 法では、以下の方程式を一次の密度行列が収束するまで繰り返す。
𝐅 X (±ω) = 𝐝X + 𝐉𝐃X (±ω) + (𝐕 XC [𝐃X + 𝐸𝐃X (±ω)] − 𝐕 XC [𝐃X])/𝐸(1)

ここで、𝐅は Fock 行列、𝐝は双極子モーメント行列、D は密度行列、ωは外場の周波数、𝐉は交換
積分、𝐕 XCは DFT 交換相関ポテンシャルであり、E は微小な外部電場である。上付きの X, Y, Z はベクトル成分、(±ω)は一次摂動項を表している。このとき、全体の Fock 行列𝐅 X (±ω)から部分系の Fock 行列𝐅 Xs (±ω)を抜き出し、変換行列 U を次式のように局在化領域で構築することによって計算コストを線形に削減することが可能となる。
Xs(±ω) =U𝑎𝑖(𝐂 0s† 𝐅Xs (±ω)𝐂0s )𝑎𝑖0sε0s𝑖 − ε𝑎 ∓ ω(2)

𝐂は係数行列、 は軌道エネルギー、上付きの 0 は無摂動項、a, i はそれぞれ非占有軌道、占有軌道である。局在化領域における変換行列𝐔 Xs (±ω)は次の式で定義される。
𝐂 Xs (±ω) = 𝐂 0s 𝐔 Xs (±ω)
この式から得られた𝐂(3)

Xs (±ω)を用いて、局在化領域に対する一次の密度行列𝐃Xs (±ω)を構築する。

𝐃Xs (±ω) = 𝐩𝐂 Xs (±ω)n𝐂0s† + 𝐩𝐂 0s n𝐂Xs† (∓ω)
Xs (±ω)を用いて全体の密度行列ここで p は分割行列、n は軌道の占有数を表す。得られた𝐃(4)

𝐃X (±ω) = ∑ 𝐃Xs (±ω)

(5)

s

を作り、
αXY(±ω) = Tr[−𝐝X 𝐃Y (±ω)](6)

から周波数依存分極率を求める。周波数依存分極率は励起エネルギーにおいて発散するので、その極から励起エネルギーを見積ることができる。

【DC-TDCPDFT の数値検証】
非局在系であるポリエン鎖 C30H32 に対して上述の DC-TDCPDFT を用いて求めた 0-6 eV における周波数依存分極率の XX、YY、ZZ 成分の平均を Fig. 1 に示す。二重結合間の二つの炭素とそれらに結合する水素からなるユニット(-CH=CH-)
Table 1. Excitation energy Eex and oscillator
を部分系とした。ただし、両端は(-CH=CH2)で strength of polyene chain C30H32 calculated byある。また、バッファ領域は左右 6 ユニットで standard TDDFT and DC-TDCPDFT.
DCある。基底関数は 6-311G**、汎関数は LC-BLYP
TDDFT
TDCPDFTを用いた。これによって求められた励起エネル
Eex
Eex
Oscillator
State
Transition
ギ ー と 、 同 じ く 汎 関 数 LC-BLYP を 用 い た
[eV]
[eV]
Strength
TDDFT によるポリエンの励起エネルギーと振 1Bu 2.96
6.08
allowed2.97
0.00
forbidden
動子強度の結果を Table 1 に示す。この 2 つを 2Ag 3.62
0.52
allowed4.26
比較すると、振動子強度が比較的大きい励起状 2Bu 4.28
0.00
forbidden
態 1Bu, 3Bu, 4Bu の励起エネルギーを誤差 0.11 2Ag 4.66
0.00
forbiddeneV 以下と精度よく記述できていることがわか 3Ag 4.88
3B
5.01
0.00
alloweduる。状態 6Bu のときの状態は極として出現して0.00
forbidden
いない。これは振動子強度が小さいため、極が 4Ag 5.42
0.17
allowed
5.54
出現する周波数の範囲が非常に狭くなった結 4Bu 5.43
0.00
allowed
果であると考えられる。また、禁制遷移に対す 1Au 5.82
0.00
allowed
る励起状態は極として出現しない。これから、 5Bu 5.83
0.10
allowed
周波数依存分極率を用いることで振動子強度 6Bu 5.86
0.00
forbiddenが大きいものに関しては励起エネルギーを正 1Bg 5.91
0.00
forbidden
しく見積れることがわかる。当日はポリエンの 5Ag 5.93
ほか、CT 励起を伴う系、局在系、溶媒系に対8000する結果もあわせて報告する。
6000
4000
Polarizability α/a.u.

[1] K. Kitaura, E. Ikeo, T. Asada, T. Nakano, and M.
Uebayashi, Chem. Phys. Lett. 313, 701 (1999).
[2] W. Yang, Phys. Rev. Lett. 66, 1438 (1991).
[3] M. Kobayashi and H. Nakai, in Linear-Scaling
Techniques in Computational Chemistry andPhysics: Methods and Applications (2011,
Springer), pp. 97-127.
[4] D. G. Fedorov and K. Kitaura, J. Chem. Phys.
123, 134103 (2005).
[5] M. Kobayashi and H. Nakai, J. Chem. Phys. 131,
114108 (2009).
[6] T. Touma, M. Kobayashi, and H.Nakai, Chem.
Phys. Lett. 485, 247 (2010)

2000
0
-2000
-4000
-6000

2.96

4.28

5.43

4

5

-8000
0

1

2

3
Frequency ω/eV

Fig.1 Frequency-denpendent
polarizability of polyene chain C30H32.

6

4P-088

FMO 法によるエストロゲン受容体の構造揺らぎを考慮した分子内相互作用解析(東海大・情教セ 1、東大・生研 2、みずほ情報総研 3)
合田（日向寺）祥子 1, 渡邊千鶴 1,2, 福澤 薫 3

Intramolecular Interaction Analysis of
Estrogen Receptor α by ab initio FMO Method(Tokai Univ.1, Univ. of Tokyo2, Mizuho Info&Res Inst.3)
Sachiko Aida-Hyugaji1, Chiduru Watanabe1,2, Kaori Fukuzawa3

１．序DNA の転写を制御するタンパク質の一種である核内受容体は、ホルモン等の低分子（リガンド）と結合することによって引き起こされるアロステリック効果によりその制御機能を示す。結合するリガンドの種類によって変化する構造は異なり、核内受容体に結合するリガンドは転写活性の亢進剤（アゴニスト）もしくは拮抗剤（アンタゴニスト）として働くことが知られ、乳癌等の疾患と深く関わっている。したがって、この構造変化の原因を解明することは、創薬において重要な知見を与える。
核内受容体の一種であるエストロゲン受容体α（ERα）において、結合するリガンドがアゴ
ニストもしくはアンタゴニストとして働くかどうかは、C 末端に存在するα-へリックス（H12）の位置（アゴニストポジション／アンタゴニストポジション）の違いが決めている。本研究では、
ERαにおける H12 の位置決定メカニズムを解明するため、H12 の各残基に働く分子内相互作用の解析を行った。
２．方法
まずは、ERαとアゴニスト（EST）およびアンタゴニスト（OHT）との複合体の X 線結晶構造を用い、フラグメント分子軌道（FMO）計算を、FMO-MP2/6-31G*レベル、真空中の条件下にて行った。FMO 計算の特徴として、フラグメント間相互作用エネルギー（IFIE）の解析を行えることが挙げられるが、本研究では H12 領域と相互作用して位置決定に寄与する ERα内のアミノ酸残基を特定するため、この IFIE 解析を行った。
なお、生体内におけるタンパク質周辺の環境は、溶媒や熱などの影響を受けて揺らいだ構造をとる。本研究では、この揺らぎを考慮すべく、MD シミュレーション（力場：AMBER99SB、GAFF）
により水溶液中における ERα－EST 複合体の揺らいだ状態を再現し、構造サンプリングを行った。これらの構造に対して、X 線結晶構造の時と同様の FMO 計算（FMO-MP2/6-31G*、真空中）

および IFIE 解析を行い、リガンド周りの水素結合ネットワークや、H12 の位置決定に寄与する
アミノ酸残基への影響を解析した。
３．結果および考察
３－１．X 線構造を用いた計算
H12 領域は負の電荷を持つため、周辺に位置する Lys と強く相互作用する傾向が見られた。
特に、ERα－EST（アゴニスト）複合体においては、Lys529 が大きな相互作用エネルギーを示し、H12 をアゴニストポジションに固定しているものと考えられる。さらに、Glu339 と Glu419
が H12 の N 末端側にある Lys531 を強く引きつけることで、H12 がアゴニストポジションに折りたたまりやすくなっていると推測される。（図１）
一方、ERα－OHT（アンタゴニスト）複合体では、リガンドサイズの影響を受けて Glu339
および Glu419 は Lys531 と相互作用できず、H12 は ERαのその他の領域から少し離れた場所に位置する。そのため Lys362 と強く相互作用することが可能となり、結果的にアンタゴニストポ
ジションに固定されるものと考えられる。（図１）３－２．揺らぎを考慮した計算
まず、12ns～50ns の MD シミュレーションでは、X 線構造から大きく構造が変化することはなかった。図２に MD のスナップショット構造を示す。
しかし、スナップショット構造の FMO 計算によりリガンド（EST）周りの水素結合ネットワークを解析したところ、X 線構造では Glu353、Arg394、His524 がリガンドを固定する働きを示していたのに対し、MD のスナップショット構造では Glu353、His524 のみが相互作用をしているという結果になった。Arg394 については、周辺に位置する Glu353 と相互作用する様子も観察され、リガンド結合については不安定な要素を含んでいる可能性が示唆された。
また、H12 と他のアミノ酸残基間の IFIE 解析では、X 線構造を用いた計算の結果同様、Lys529が特に強く H12 と相互作用しており、H12 をアゴニストポジションへ固定するメカニズムは生体内でも安定であるものと考えられる。

Lys362(アンタゴニスト)

Helix 12

EST

Lys529(アゴニスト)

図１：H12 と相互作用する Lys

図２：MD のスナップショット構造

4P089CPMD を用いた梯子型 π 共役化合物の発光に関する理論的研究(京大院工 1，京大福セ 2) 成田 真弘 1，佐藤 啓文 1，中尾 嘉秀 2

Theoretical Study of Luminescence of Ladder-type π-conjugatedCompounds by the CPMD method(Kyoto Univ.1) Masahiro Narita1, Hirofumi Sato1, Yoshihide Nakao1【序】
多環式骨格を含む Scheme 1 に示すような梯子型 π 共役化合物は平坦で剛直な分子構造を持ち、安定で秩序立った結晶構造をとることが知られている。また、π 共役系に様々なヘテロ元素を導入することで電子状態を、側鎖 R を変えることで結晶構
Scheme 1造を制御できるため、新しい有機エレクトロニクス材料として非常に注目されている[1]。しかしながら、固体中の π 共役化合物の発光に関する理論計算はほとんど行われておらず、設計と応用の面からも電子状態の解明が求められている。本研究では、梯子型 π 共役化合物についてヘテロ元素の寄与を明らかにし、基底状態と励起状態の安定構造を求めることで吸収・発光エネルギーを明らかにした。また、Car-Parrinello Molecular Dynamics(CPMD)[2]計算を用いて各状態における安定構造と励起状態の構造緩和する過程を比較した。
【計算方法】
気相中の計算には基底関数として cc-pVDZ を用いた。構造最適化は B3LYP 法及び TD-B3LYP法で行い、励起エネルギーは TD-BLYP 法と TD-B3LYP 法、GMCPT 法で求めた。計算パッケージは GMCPT 法のみ GAMESS を用い、他は全て Gaussian 09 で行った。固体中の構造最適化は CPMD
パッケージで平面波基底(cutoff = 90 Ryd.)、擬ポテンシャル及び周期的境界条件を用い、BLYP 法で行った。ただし、固体中の励起エネルギーは気相中と同様に Gaussian 09 で TD-B3LYP 法によって評価した。CPMD 計算は 1step を 4 a.u. ( = 0.097 fs ) として、仮想電子質量 400 a.u.、各状態の安定構造から温度 50±2 K と 300±2 K でスケーリングし、9.7 ps 平衡化した後に 9.7 ps 分のトラジ
ェクトリを用いた。構造緩和する過程は 0.3 ps 分のトラジェクトリ 100 本から解析を行った。
【結果と考察】R = H の分子の気相中における吸収・発光エネルギーについて、TD-BLYP 法、TD-B3LYP 法、
GMCPT 法で求めた結果を Table 1 に示した。この吸収・発光には HOMO-LUMO 間の遷移が大きく寄与しているため、 状態の安定構造はキTable 1. R = H の気相中の吸収・発光エネルギー(eV)

ノイド構造となり、発光エネルギーは約 0.4

TD-

TD-

BLYP

B3LYP

Absorption

2.60

2.96

Emission

2.26

2.53

eV と大きくストークスシフトする。実験では R = Pr の THF 中の吸収エネルギーは 2.94
eV であり、TD-B3LYP 法が比較的よく一致した。TD-BLYP 法では吸収エネルギーを小[1]

a. R = Pr の THF 中の実験値

GMCPT

Exp.[a]

3.11

2.94

2.71

2.81

さく見積もり、計算コストを抑えるため励起配置に制限を課した GMCPT 法では吸収エネルギーを大きく見積もった。従って、吸収・発光エネルギーの評価には TD-B3LYP 法を用いることにした。
ケイ素原子が π 共役系に与える影響を調べたところ、LUMO で σ*(Si-H)の反結合性相互作用がπ 共役系に寄与し、吸収・発光エネルギーを約 0.2 eV 低下させることが分かった。また、ケイ素同士の位置が近い方が吸収・発光エネルギー共に小さくなることが明らかになった。
R = Pr では結晶構造が三斜晶であり、CPMD パッケージでは 状態の構造最適化ができない上に振動子強度も得られない。また、TD-BLYP 法による計算しか行えず励起エネルギーを過小評価する。そこで、固体中の 状態の構造に 状態の構造を代用し、気相中と同様に Gaussian 09 で
TD-B3LYP 法を用いて励起エネルギーを評価することにした。 構造の妥当性を確認するため、気相中における R = H について、 状態の安定構造が Gaussian 09 と CPMD パッケージでほぼ同じであること、 状態及び 状態の構造が大きく変わらず、吸収・発光エネルギーも近い値をとることを確認した。また、固体中の分子が隣接分子から受ける寄与について調べたところ、励起エ
ネルギーはほとんど変化しないことを確認した。
固体中の吸収・発光エネルギーを Table 2 にまとめた。固体中の励起状態は気相中と同様にキノイド型を安定構造とし、固体中の発光エネルギーは約 0.5 eV と大きくストークスシフトする。
Scheme 1 に示した C1-C2 結合は基底状態で二重結合性を、励起状態では単結合性を示して特に大きく変化する。そこで、気相中の分子について CPMD 計算を行い、安定構造における C1-C2 結合長の分布を調べた結果を Fig. 1 に示した。結合長の平均値は気相中における基底状態と励起状態の安定構造の結合長 1.38 Å と 1.44 Å にそれぞれよく一致し、温度の上昇に従って結合長の分布が広がる。励起状態における構造緩和過程について 100 本のトラジェクトリを解析したところ、50K
では励起状態における最安定の結合長に 7 fs 程度で緩和するという結果が得られた。また、C1-C2結合の振動モードは余剰な熱を 20 fs 程度で失う
S0 (50K)

ことが分かった。

S0 (300K)

緩和過程の詳細な比較については当日報告す

T1 (50K)

存在確率る。

T1 (300K)

Table 2. 固体中の吸収・発光エネルギー (eV)In solid

Exp. in solid

Absorption

2.90

―

Emission

2.41

2.61

1.3

1.4

1.5

1.6

結合長(Å)
Fig.1. 結合長の分布【文献】[1] H. Yamada, C. Xu, A. Fukazawa, A. Wakamiya, S. Yamaguchi, Macromol. Chem. Phys., 210, 904(2009).
[2] R. Car and M. Parrinello, Phys. Rev. Lett., 55, 2471 (1985).

4P-090
Elongation 法によるヘテロナノチューブの電子状態に関する理論的研究（九大院・総理工 1、 SCNU2、JST-CREST3）
筒井亮太 1、Liu Kai1、折本裕一 1、Feng Long Gu2,3、青木百合子 1,3

A theoretical study on electronic structures of heteronanotubes
by elongation method(Kyushu Univ1., SCNU2, JST-CREST3)
Ryota Tsutsui1, Liu Kai1, Yuuichi Orimoto1, Feng Long Gu2,3, Yuriko Aoki1,3【序】 カーボンナノチューブ(Carbon
nanotubes:CNTs)をヘテロ元素で置換した構造を持つ
ヘテロナノチューブは、CNTs とは違った物性を持つ。
特に隣合う炭素原子を窒素とホウ素で交互に置換した構造を持つ BN ナノチューブ(Boron Nitride
nanotubes:BNNTs, Figure 1)は、CNTs に比べて高い
熱的安定性や化学的安定性を持ち、またその電気的特性は極めて異なっている。CNTs が金属及び半導体的な性質を示すのに対して、BNNTs は基本的に約 5.5eV のバ

Fig. 1 (4,4) BNNT(armchair 型,上)と(6,0) BNNT(zigzag 型,下)の構造

ンドギャップを持つ絶縁体である[1]。このような特性の解析のために、量子化学的なアプローチが行われてきたが、軸方向長さがマイクロメートルオーダーにまで達する
巨大系の電子状態解析は、従来の ab initio
な方法では困難である。当研究室ではこのような巨大系に対する、高精度かつ高効率な計算手法として、ab initio 理論的重合法-Elongation 法-を開発し、改良を重ねてきた[2]。本研究では、Elongation 法を様々な構造の BNNTs のエネルギー計算に対して適用し、その電子状態の解析を試みた。
【理論】 Elongation 法は高分子の重合反応のアイデアを元にしている。Figure 2 にElongation 法の手続きの概略図を示す。まず、任意の大きさの出発クラスターに対して Hartree -Fock 方程式を解くことで、正Fig. 2 Elongation 法における手続きの概略図。
CMO は正準分子軌道、RLMO は領域局在化軌道

準分子軌道(CMOs)を得る。次に得られた CMOs に対してユニタリー変換を行い、予め定めた
Frozen 領域と Active 領域に軌道を局在化させて領域局在化軌道(RLMOs)を得る。Attacking
monomer となるフラグメントは Frozen 領域から十分に離れているので、Active 領域とフラグメントの相互作用のみを考えても全エネルギーに影響を与えない。すなわち、Fock 行列の次元を落とした計算が可能になり、計算が高速化する。伸長反応を進めていっても相互作用領域の大きさは変わらないので、サブ O(N)スケールでの計算が可能となる。さらなる高速化の手法として、当グループで開発された Cutoff-Elongation 法がある[3]。この方法では、フラグメントから一定以上
距離の離れた Frozen 領域の軌道を Cutoff することで高速化を実現している。
【結果】

Table 1 (4,4)BNNT における
Elongation 法と従来法のエネルギー計算の比較(4,4)と(6,0)のカイラリティを持つ 2 種の BNNT について計算を行った。電子状態計算は

Number of Number of
units
atoms

RHF/STO-3G レベルで行い、1 ユ

6
7
8
9
10
11
12
13
14
15
16
17
18
19

ニットはそれぞれ B8N8、B12N12とした。出発クラスターを N=6 とし、ユニットを１つずつ増やして
いき N=19 まで伸長させた。N=13となった時点で軌道のカットオフを行っている。Elongation 法の計
算精度の検討のために、各ステッ
プにおいて同様の設定で従来法ににおけるエネルギー計算の結果を示した。EELG は Elongation 法におけるエネルギー、ECONV は従来法におけるエネルギーで、E/atom を(EELG – ECONV)/原子数と定義する
ことで、計算精度の指標とした。
(4,4)及び(6,0)の両構造において
E/atom の値は 10-8 から 10-10 のオーダーになっており、非常に高い精度で計算を実行できている事がECONV-3765.7529734821-4391.9389523565-5018.1249306711-5644.3109086737-6270.4968864941-6896.6828642046-7522.8688418462-8149.0548194427-8775.2407970087-9401.4267745561-10027.6127520897-10653.7987296138-11279.9847071314-11906.1706846441

E/atom(a.u.)
0.000E+00
2.312E-10
4.604E-10
6.863E-10
8.985E-10
1.095E-09
1.268E-09
2.919E-09
3.045E-09
3.080E-09
3.111E-09
3.145E-09
3.175E-09
3.203E-09*[4]

Table 2 (6,0)BNNT における
Elongation 法と従来法のエネルギー計算の比較

おける計算を行った。
Table 1 及び 2 にそれぞれの構造

112
128
144
160
176
192
208
224
240
256
272
288
304
320

Energy (a.u.)
EELG-3765.7529734821-4391.9389523269-5018.1249306048-5644.3109085639-6270.4968863360-6896.6828639944-7522.8688415825-8149.0548187888-8775.2407962779-9401.4267737676-10027.6127512435-10653.7987287080-11279.9847061663-11906.1706836192

Number of Number of
units
atoms
6
156
7
180
8
204
9
228
10
252
11
276
12
300
13
324
14
348
15
372
16
396
17
420
18
444
19
468

Energy (a.u.)
EELG
ECONV
-5641.7295638111-5641.7295638111-6580.9347484942-6580.9347487575-7520.1399627586-7520.1399636792-8459.3451953119-8459.3451972159-9398.5504396528-9398.5504427921-10337.7556918767 -10337.7556964400-11276.9609495434 -11276.9609556746-12216.1662113133 -12216.1662188861-13155.3714756982 -13155.3714850069-14094.5767421727 -14094.5767532993-15033.7820102533 -15033.7820232522-15972.9872795850 -15972.9872944943-16912.1925498986 -16912.1925667532-17851.3978210053 -17851.3978398310

E/atom(a.u.)
0.000E+00
1.463E-09
4.513E-09
8.351E-09
1.246E-08
1.653E-08
2.044E-08
2.337E-08
2.675E-08
2.991E-08
3.283E-08
3.550E-08
3.796E-08
4.023E-08*[4]

分かる。このことから、BNNTs の電子状態解析において、Elongation 法が信頼出来る計算手法であると言える。発表当日はより大きな基底関数での検討や、電場に対するエネルギーの応答を見る Finite-Field 法を当手法に応用した Elongation-FF 法により BNNTs の NLO 特性に関する解析を行った結果を示す予定である。
【参考文献】[1] D. Golberg, Y. Bando, C. Tang, and C. Zhi, Adv. Mater. 19, 2413-2432 (2007)[2] A. Imamura, Y. Aoki, and K. Maekawa, J. Chem. Phys. 95, 5419-5431 (1991).
[3] J. Korchowiec, F. L. Gu, A. Imamura, B. Kirtman, and Y. Aoki, Int. J. Quantum Chem. 102,785-794 (2005)[4] W. Chen, G. Yu, F. L. Gu, and Y. Aoki, J. Phys. Chem. C 113, 8447-8454 (2009).

4P091
量子輸送現象におけるシアニドおよびイソシアニドの	 アンカーとしての影響	 （九大先導研）古賀 惇也，辻 雄太，吉澤 一成

Orbital Control of Single-Molecule Conductance Perturbed by π-accepting Anchor Groups:
Cyanide and Isocyanide(IMCE Kyushu Univ.)

Junya Koga, Yuta Tsuji, and Kazunari Yoshizawa【序】単一分子デバイスは近年限界が指摘されつつあるシリコンベースの半導体デバイスにとって代わる次世代のデバイスとして注目を集めている[1]。単一分子デバイスにおいて非常に重要な要素としてアンカー部位がある。アンカー部位は、単一分子を電極と接合させる役割を担い、この部位が電子輸送物性に与える影響は大きい。代表的なアンカー部位として、金と強く相互作用するチオール基 (-SH) が広く採用されてきた。しかし、吸着構造の変動性に関する問題などが指摘され、他の官能基を用いた研究が近年のトレンドとなっている[2]。本研究では、強力な電子吸引基であるシアニド (-CN) およびイソシアニド (-NC) をアンカー部位として用いた場合の単一分子の電気伝導度特性をヒュッケル法および密度汎関数法レベルの非平衡グリーン関数法を用いて理論的に検討し、これらのアンカー部位が電子輸送物性に与える影響を分子軌道の観点から明らかにすることを目的としている。

【計算方法】近年，フロンティア軌道の位相と振幅から

N

C

C

N

単一分子の電子輸送物性を定性的に予測する規則が提案されている[3]。効果的な電子輸送が起きるためには以下の 2 つの条件を満足する必要がある。(1) HOMO とLUMO の振幅が大きい原子を電極に接続しなければならない，(2) 電極に接続される 2 原子上の分子軌道係数

N

C

N

N

C

p-BDCN

p-BDNC

C

C

N

C

N

N

C

の積の符号が HOMO と LUMO で異なっていなければならない。アンカー部位としてチオール基を用いた場合においても同様の予測が有効であることが報告されていm-BDCNm-BDNC

図 図１!計算モデル
1. 計算モデルる[4]。
	
  本研究では、図 1 に示した 4 つの分子に対し解析を行った。まず、これらの分子のヒュッケル
分子軌道に上の規則を適用し、それぞれの電子輸送について定性的予測を行った。続いて、ヒュッケル法と非平衡 Green 関数法を組み合わせた NEGF-HMO 法を用いた電気伝導の計算を行い、定性的予測との比較を行った。定量的な解析として、密度汎関数理論と非平衡 Green 関数法を組み合わせた NEGF-DFT 法を用いた計算を行い、定性的予測との比較を行った。

【結果と考察】	  図 2 に BDCN のパラ体とメタ体のヒュッケル法で計算したフロンティア軌道を示す。上述の規則をこれらの軌道に適用するとパラ体では HOMO と LUMO のコンダク
タンスへの寄与は互いに強め合うもので、メタ体では打ち消し合うものであることが分かる。従って、パラ体はメタ体に比べ

図 2. BDCN のフロンティア分子軌道高いコンダクタンスを持つと予測できた。これは、BDNC の場合でも同様であった。
	  図 3 に BDCN のパラ体とメタ体に対し、NEGF-HMO 法および NEGF-DFT 法で計算した電子の透過確率を示す。図 3(a)では、フェルミエネルギーの位置 (E = 0) で、メタ体の透過確率が 0 へと落ち込んでいることが分かる。また、図 3(b)においても、フェルミエネルギー近傍にて透過確率はパラ体の方がメタ体に比べ、大きい値を示した。同様の結果が、BDNC の場合においても得られた。これらの結果はヒュッケル分子軌道を用いた定性的予測とよく一致する。フロンティア
軌道の位相と振幅を用いた定性的なコンダクタンスの予測が、シアニドおよびイソシアニドをア
ンカー部位として用いた場合においても有効であることが明らかとなった[5]。

図 3. (a) NEGF-HMO 法および (b) NEGF-DFT 法により求めた電子の透過確率[1] Chen, F.; Tao, N. J. Acc. Chem. Res. 2009, 42, 429.
[2] Li, Z.; Borguet, E. J. Am. Chem. Soc. 2012, 134, 63.
[3] Yoshizawa, K.; Tada, T.; Staykov, A. J. Am. Chem. Soc. 2008, 130, 9406.
[4] Tsuji, Y.; Staykov, A.; Yoshizawa, K. J. Am. Chem. Soc. 2011, 133, 5955.
[5] Koga, J.; Tsuji, Y.; Yoshizawa, K. J. Phys. Chem. C, under revision.

4P092
4P092
Quantum chemical DMRG study of biologicalbiological dinuclear metal
metal
2+
9
complexes: ∆ desaturase and {{[[Cu(NH3)3]2O2} complex(Institute for Molecular Science1) Jakub Chalupský1, Yuki Kurashige1, Takeshi Yanai1

We present quantum chemical study of two biologically relevant dinuclear transition
metal complexes based on DMRG-CASSCF and DMRG-CASPT2 calculations. Thesemethods, allowing us to use sufficiently large active spaces for dinuclear complexes, canprovide useful information about electronic structure of low-lying states as well as theenergetics of studied processes. We address two problems in this poster, barrier hight of
one of the reactions in catalytic cycle of ∆9 desaturase and energetics of isomerization of{[Cu(NH3)3]2O2}2+ complex.
Stearoyl-ACP ∆9 desaturase is the non-heme iron-containing enzyme involved in fattyacid metabolism and its regulation in plants, which catalyzes the conversion of thesaturated stearic-acid chain into the unsaturated oleic-acid chain. It undergoes multistep catalytic cycle, whose mechanism is under investigation by the experiments as wellas theoretical DFT-based QM/MM calculations. Although several alternatives of its
mechanism have been already revealed by theory, reliable prediction of the barrierheight for one of the key steps – abstraction of the first hydrogen atom from thesaturated stearic-acid chain – seems to be very difficult to make using DFT or
conventional CASSCF/CASPT2 approaches. We thus provide a study of height of thisbarrier based on DMRG-CASSCF/CASPT2 methodology, which allows for qualitatively
correct description of the electronic structure of the active site. We have found that thistype of calculations is practically feasible even for such a large system (model of active
site involves almost 70 atoms), and it provides very important results. Our DMRGCASPT2 calculations have shown that from the two alternatives of the mechanism
studied so far, one (so called “protonated” version of the mechanism, geometry change

during the studied step is shown in the figure) is significantly more favored than theother (water-containing version). Although a lot of work still needs to be done for otheralternatives of the mechanism of catalytic cycle of ∆9 desaturase, already these firstobtained results seems to be very promising, and we expect this project to be a highly
important application of DMRG approach to biologically relevant chemical problems.
Another interesting problem is the energy of isomerization between bis(µ-oxo) andperoxo forms of {[Cu(NH3)3]2O2}2+ complex, which are mimicking structural motives invarious copper-containing enzymes, for example tyrosinase. Accurate predictions of therelative energies of these structures have been shown to be extremely difficult to makeby nowadays theoretical methods – they seem to require either RASPT2 with the size ofactive space practically at its limit, or CCSD(T) including corrections formulticonfigurational character of the reference wave function. We thus consider the
problem of energetics of isomerization of {[Cu(NH3)3]2O2}2+ complex (in the figure) to be
an important case for testing the applicability of DMRG approach on systems, which

require including many orbitals on one center in the active space (Cu 3d and 4d shell –10 orbitals on a single atom), because it is expected from the principles of DMRG that
this type of systems might be problematic to describe correctly. This study thus should
provide important information about way of dealing with such systems and possible
limitations of DMRG-CASSCF/CASPT2 approach, which could be of high importance tothe theoretical chemists in this field.

4P093

DNA 伸長反応についての ONIOM 計算（筑波大院・化）廣田 恵以、守橋 健二

ONIOM calculation on DNA elongation processes(University of Tsukuba) Kei HIROTA, Kenji MORIHASHI【序論】
DNA ポリメラーゼによる一本鎖 DNA 鎖の伸長過程で重要となる核酸塩基対の相互作
用は興味深い計算対象である。核酸塩基対によるπ－πスタッキングの存在は、伸長反応を行う上で大きく貢献すると考えた。そこで本発表では、ONIOM 法を用いた DFT 計算を行い、塩基間相互作用が大規模系での伸長反応にどのように寄与するかを明らかにした。DNA 鎖伸長反応については、ポリメラーゼの二つの Mg イオンの関与によるリンの五配位中間体を介した反応経路を採用した。核酸塩基対は cytosine-cytosine 対を計算し、この核酸塩基対を QM 部位に含むか否かで結果がどのように変化するかを比較、検討した。
【計算対象】DNA 鎖伸長反応には Polymeraseβ(Polβ)が関与する。そこで、X 線結晶構造[1](PDB file: 2fmp)を基にして活性部位をモデル化したもの（図 1）を考え、これらを計算対象とした。
R1 は本来 DNA が結合しており、今回のモデル化では糖五員環と核酸塩基対を1 つだけ付加させ、残りは省略したものとなる。R2 も同様に糖五員環と核酸塩基対を付加して計算をしており、この R1

図 1. Polβ複合体モデルの構造。R1、R2 には核酸塩基が結合している。ここから、primer の酸素原子が Pαへと求核攻撃を行い、反応が進行してと R2 の核酸塩基対間でπ－πスタッキングが生じている。
本反応の酵素となる Polβは 2 つの Mg2+イオンを含んでおり、これが六配位構造を取ることで求核攻撃をより起こりやすくしている。
【計算方法】
本反応における反応経路を DFT 計算により求め、本研究では汎関数系としてπ－π
スタッキングを評価できる M06-2X[2]を用いた。作成したモデルの総電荷は-3 として計算し、計算手法は M06-2X/6-31++G(d,p)//M06-2X/6-31G(d,p)で行った。また、比較対象
として核酸塩基対を MM 部位とした計算も行い、計算手法は M06-2X/6-31++G(d,p)//
ONIOM (B3LYP/6-31G(d):UFF)で行った。
本反応は primer の H が Pαの O に移り、Pγの O へと移動していく。中間体では Pαを中心とした五配位中間体となっている。また遷移状態を経由して、終状態では primer
からピロリン酸が解離する。各計算手法において得られた相対エネルギーを図 2 に示した。
【結果・考察】始状態を基準とした計算の結果、中間体の構造が変化した。それに伴い、相対エネル
ギーの値も大きく変化をした。構造の変化は図 2 に示す通りである。

中間体(ΔE＝1.90)

始状態 (ΔE＝0.00)

終状態(ΔE＝5.96)
さ差中間体(ΔE＝15.63)
図 2. 計算手法による反応経路の比較 (kcal/mol)
黒丸で囲んだプロトンが転移することで反応が進行する。

このように中間体の構造が大きく異なり、活性化エネルギーについても変化が見られた。
このことから、π－πスタッキングの存在は計算手法によって大きく左右することが確かめられた。また大規模系に着手する際についても、π－πスタッキングによる活性障壁の減少が無視できないものであるため、QM 部位の計算範囲についても大きく変化させる必要が出てくるのではないかと考えられる。
[1] V. K. Batra et al., Structure, 14, 757 (2006).
[2] Y. Zhao and D. G. Truhlar, Theor. Chem. Acc, 120, 215 (2008).

4P094RISM 法を用いた N-acetylproline amide の円偏光二色性スペクトルに関する理論的研究（神戸大院・シス情 1，理研 AICS2） 北川 裕也 1，秋永 宜伸 2，天能 精一郎 1Theoretical study of the circular dichroism spectra of N-acetylproline amideusing RISM method（Kobe Univ.1，RIKEN AICS2） Yuya Kitagawa 1，Yoshinobu Akinaga 2，Seiichiro Ten-no 1【序】
円偏光二色性(CD)スペクトルは分子の左右円偏光に対する吸光度の差に相当する物理量であり、分子の三次元構造の変化を鋭敏に反映する。このためキラル分子の同定やタンパク質の折り畳みなど、分子の構造解析に広く用いられている。また CD スペクトルは分子の構造だけでなく溶媒の種類など周辺環境にも強く依存するため、理論計算においては周辺環境も考慮する事が重要である。今回、ab initio 法に RISM 法を組み合わせることで、水溶液中の分子の CD スペクトルの計算を行なった。そして気相中及び水溶液中の分子の構造と溶媒効果について調べ、CD スペク
トルにおける RISM 法の有効性を調べた。
対象分子は N-acetylproline amide(NAPA)である。NAPA はアミノ酸のプロリンの誘導体であり、もっとも単純なペプチドの周辺環境への依存性を調べるためのモデルの一つである。NAPA
は Fig.1 の構造をしており、2 つの二面角γ、ψにより 2 つのアミド基の配向が決まり、cis 及びtrans 型の C7、PII、310 と呼ばれる配座異性体が生じる。これら

H2
Cはさらに二面角φの違いからそれぞれにアキシャル(Ax)/エカト
H2C

リアル(Eq)配座が存在しうる。NAPA は CD スペクトルなどの分φN

光実験から非極性溶媒中と水溶液中ではコンフォメーションが異なると考えられており、水溶液中の構造については現在も議論が行われている。これまでの NAPA に関する研究では Eq 配座のみ

C H2

H3C

CH

γ ψ
O

NH 2

O

を扱ったものであったが、今回は Ax 配座を含めて CD スペクト

Fig.1 N-acetylproline amide

ルへの寄与を調べた。

(trans PII 配座)【方法】
複数の NAPA の初期構造から、気相中及び水溶液中(PCM)にて MP2/cc-pVDZ で構造最適化を行なった。構造最適化により得られた各配座について、CC2/aug-cc-pVDZ で気相中、水溶液中の励起状態及び回転強度(Rn)の計算を行なった。Rn の計算にはゲージ原点に依存しない速度表現を
用いた。

RnV =

1
2ωn

Ψ 0 p̂ Ψ n ⋅ Ψ n ˆl Ψ 0水溶液中における溶媒の影響は RISM-PT1 法[1]により計算した。パラメータは OPLSAA のTIP3P を使用した。

【結果・考察】
気相中において C7-Ax/Eq、310-Ax/Eq など 7 つの安定な配座が得られた。このうち最も安定な C7 配座について CD スペクトルの計算を行なった。計算および実験スペクトルを Fig2 に示す。C7-Ax、C7-Eq のCD スペクトルは共に 230 nm、190 nm 付近に負のピークが現れよく似た特徴を示した。しかしながらEq 配座では 205 nm 付近に正のピークが見られ、実験スペクトルと一致しない。一方、Ax 配座のスペク

Fig.2 気相中の CD スペクトル

トルは実験スペクトルとほぼ一致した。安定構造のエ
ネルギーを比較すると Ax 配座の方が 1 kcal/mol 程度低いことから気相中では C7-Ax 配座がより支配的に存在し CD スペクトルへ寄与していると考えられる。
プロトン性極性溶媒中における NAPA の構造については未だ議論が行われている。水溶液中での構造最適化により C7-Ax/Eq、310-Ax/Eq、PII-Ax/Eq など 10 の配座が得られた。これらのうち
trans 体の C7、PII、310 の Ax 及び Eq 配座について RISM-PT1(TIP3P)で溶媒効果を取り込んでCD スペクトルの計算を行なった。得られた CD スペクトルおよび実験スペクトルを Fig3 に示す。
実験による NAPA の水溶液中の CD スペクトルは気相中などの非極性媒質中に比べ強度が弱い負のピークが見られるのが特徴である。気相中とは異なり水溶液中での計算では単一の配座で実験
スペクトルに一致するものは見られない。これについては水溶液中では配座間の自由エネルギー差が小さくなり NAPA が有意に複数のコンフォメーションを取っていることが示唆される。
NAPA の気相中と水溶液中の CD スペクトルの違いが何に起因するのかを明らかにするため、水溶液中で最適化した PII-Eq の構造を用いて気相中の CD スペクトルを計算した。Fig.4 に計算結果を示す。気相中ではカルボニル基の nπ*遷移に基づくピークの赤方シフトや他の遷移についても変化が見られた。その結果、気相中で計算したものでは 235 nm、190 nm 付近に負のピーク、
210 nm 付近に正のピークが現れた。このことから NAPA の気相中と水溶液中の CD スペクトルの違いは、コンフォメーションの変化と共に溶媒との相互作用による影響に起因していることが分かる。詳細については当日に発表を行う。

Fig.3 水溶液中の CD スペクトル

Fig.4 PII 配座における気相中と水溶液中のCD スペクトルの比較

References[1] S. Ten-no, J. Jung, H. Chuman, and Y. Kawashima, Mol. Phys., 108, 327 (2010).
[2] Y. Kitagawa, Y. Akinaga, Y. Kawashima, J. Jung, S. Ten-no, Chem. Phys., 401, 95 (2012).
[3] V. Madison and J. Schellman, Biopolymers, 9, 511 (1970).

4P-095

中性条件に基づく新規静電ポテンシャルの開発（近畿大学、先端技術総合研究所 高圧力タンパク質研究センター）米澤 康滋

A long range electrostatic Potential relied on charge
neutrality(High Pressure Protein Research Center, Institute of Advanced Technology, Kinki
University)
Yasushige Yonezawa【序】分子シミュレーションにおいて長距離静電クーロン相互作用の寄与は大変大きくシミュレーション結果を左右することが知られている。すなわち分子シミュレーションおけるクーロン相互作用の取り扱いは極めて重要である。クーロン相互作用の有効作用領域は無限大であるが実際の分子シミュレーションで無限遠からの寄与を取り扱うことは不可能であり、その効果の適切な評価に関しては現在まで様々な手法が開発されて来た。比較的小さな孤立系ではすべての対象原子及び分子を考慮することが可能であるが、大きな孤立系では cutoff 法やマルチポール展開法等が使われている。また周期系では無限に続く周期系の影響を考慮し得る Ewld 法やその高速フーリエ変換を用いた高速計算版である Particle Mesh Ewald（PME）
法が標準法として頻用されている。
近年、計算機性能の向上に伴って蛋白質などの大規模分子シミュレーションが行われるよ
うになって来た。これらの大規模分子シミュレーションは顕な溶媒分子を含むことが多くそれらの系の原子数は数百万以上に及ぶこともある。このような大規模系の周期的境界条件下での計算は通常、並列計算機を用いて計算される。ところが PME 法で用いられている高速
フーリエ変換法は、大規模な並列計算には対応が難しいという困難を抱えている。また cutoff
方は対象原子に相互作用する原子を球対称な有限範囲内に制限して計算する方法で大規模並列計算に大変適している手法であるが、この方法を単純に分子シミュレーションの長距離クーロン相互作用に適応すると深刻なアーティファクトを生じて信頼性のある分子シミュレー
ションを遂行できないことが良く知られている。
Wolf とその共同研究者は、1999 年に新規なクーロン相互作用計算法を提案した。彼らの方法は、cutoff 法計算時に cutoff 半径内の電荷がゼロとなる時にそのエラーが大変小さく抑制される観測事実を基に考案された。Wolf 等は、中心となる対象原子に関して cutoff 半径内でクーロン相互作用する電荷の延長線上の cutoff 球上に反対電荷を持つカウンター点電荷を置いて cutoff 半径内の電荷を見かけ上無くす事でこの電荷中性条件を満たす方法を考案した。

さらに彼らは電荷中性条件だけでは消せない cutoff 長に依存するポテンシャル振動を、補誤差関数でポテンシャルダンピングすることで抑制した。彼らのポテンシャル（以下、Wolf ポ
テンシャル）は孤立系や周期系に適用可能で、生体高分子系を含むクーロン相互作用する様々な系で Ewald 法や PME 法の結果を安定して再現できることが検証されている。
本発表では、Wolf の中性電荷条件のみに基づく新規なクーロン相互作用計算方法について報告する。この新規クーロン相互作用法は Wolf ポテンシャルで使われた補誤差関数によるポテンシャルダンピングを必要としないので、従来の Wolf ポテンシャルによる計算法よりも高速な計算が可能である。
【理論・計算実験】
Wolf ポテンシャルと本報告の新規ポテンシャルを以下に示す。

 erfc( rij ) erfc( rij )   erfc  RC    N 2
 lim 

EWolf    qi q j   
  qirij  RC
rij
rij
2 RC
  i 1i 1 j ( i ) 
r R
N

ij

C

 1 erf ( Rc ) 1
 N 2 erf ( Rc )1 N
ETotal    qi q j  
 2 ( Rc  rij )    qi
2 i 1 j i
r
R
Rc
Rc
c
 ij
 i 1rij  Rc

EWolf は Wolf ポテンシャルの、ETotal は本報告の新規ポテンシャルで計算される全エネルギーである。N は電荷を持つ全原子数を示す。erfc は補誤差関数、erf は誤差関数である。Rcは cutoff 長を表す、α はポテンシャルのパラメーターで Wolf 法に於いてはポテンシャルのダ
ンピングファクターを示し、本報告の新規ポテンシャルではカウンター電荷の分布度合いを表している。本報告ポテンシャルの誤差関数項は、Wolf の電荷中性条件の物理的意味を深く考慮することで Wolf 等が使用したカウンター補正点電荷を、広がりのあるガウス型電荷に変更した結果得られたものである。本報告の新規ポテンシャルでは、計算負荷の高い誤差関数（補誤差関数も同様）を計算の初めに一回のみ計算して保持しておくだけで良い、一方、Wolf
ポテンシャルではポテンシャル計算毎にその負荷の高い補誤差関数の計算を繰り返す必要がある。
【結果と考察】
本新規ポテンシャルを様々な系に適用して、PME 法と比較することでその評価を行った。
その結果、本新規ポテンシャルは安定な分子シミュレーションを実現しかつ、PME 法の結果を大変よく再現できることを示すことができた。計算結果の詳細と考察については当日詳し
く報告を行う。

4P-096
CD 法に基づくタンパク質全電子カノニカル分子軌道計算プログラムの開発(東大生研)平野 敏行，佐藤 文俊

Development of canonical molecular orbital calculation program for proteinsusing Cholesky decomposition(Institute of Industrial Science, the University of Tokyo)
Toshiyuki HIRANO, Fumitoshi SATO
【概要】
分散メモリ型並列計算機において、大規模タンパク質の全電子カノニカル分子軌道計算を行う
ためには、計算に必要な大規模行列を各計算ノードに分散保持する必要がある[1]。Kohn-Sham
行列におけるクーロン項ならびに Fock の交換項の行列要素を求める際に、計算律速ルーチンの一つである 4 中心 2 電子積分が必要であるとともに、疎とはみなせない密度行列要素を参照しなければならない。したがって、大規模分子の全電子カノニカル分子軌道計算では、高速かつ高度に並列化された 4 中心 2 電子積分計算に加え、分散保持した大規模な密度行列に対する並列性能の良いアクセスが求められる。
本研究では、次世代スーパーコンピュータ「京」をはじめとする、分散メモリ型並列計算機に適したカノニカル分子軌道計算プログラムの開発を目的として、コレスキー分解[2]を用いた並列計算法を開発・実装した。巨大メモリを必要とするコレスキー分解法であるが、その行列を縮約し、また大規模行列を分散保持することにより、その欠点を克服することができた。
【コレスキー分解による大規模電子状態計算】
原子軌道のペアの次元を持つ、4 中心 2 電子積分を行列要素とする行列 は正値対称行列であり、コレスキー分解によりランク

ここでの行列 に分解できる。

は原子軌道のインデックスをあらわす。このときコレスキー分解による計算誤差は

以下のように求められ、任意の閾値 に対して計算精度を設定することができる。

ここで のサイズはとなり、巨大な行列になる。そこで本研究では、CDAM 法[3]とピボット CD 法[4]を組み合わせることにより、行列サイズを縮小した。すなわち CDAM 法では、以下の条件式
で生き残った原子軌道ペア()から構成される行列 に対してコレスキー分解を行う。これはSchwartz の不等式でのスクリーニング[5]に相当する。コレスキー分解は、ピボット CD 法により、なるべくランクが少なくなるようにコレスキーベクトルを求める。計算コストの高い 2

電子 4 中心積分は、コレスキー分解処理に必要なタイミング(on-the-fly)で求めることができるた
め、余計な積分計算が発生しない。なお、計算途中で必要なコレスキーベクトルは、行ベクトルとして各プロセスに分散して保持している。これにより、計算可能なサイズが、各計算ノードではなく、計算機全体に搭載されたメモリ容量に依存することになるとともに、ロードバランスの均等化を図っている。
コレスキー分解法において、クーロン項 ならびに Fock の交換項 は以下のように求められる[3]。

SCF 繰り返し計算におけるクーロン項 ならびに Fock の交換項 の計算において、シンプルな行
列演算のみで求めることができ、並列計算機に適した、ScaLAPACK[6]などの線形演算ライブラ
リを利用することができる。均等なロードバランスを取ることが難しい、大域密度行列のアクセスと分子積分計算が必要ないことも長所の一つである。
本手法によるコレスキー分解法では、計算精度を左右する閾値は、 と の 2 つ存在する。適切な閾値を決定するために、いくつかの分子において、これまでの update 法を用いた SCF 計算と本手法を比較した。その結果、分子積分のカットオフ値を度の計算精度を保つためには、、に設定した計算結果と同程が妥当であることがわかった。

当日は、上記コレスキー分解を用いた電子状態計算法の各アルゴリズム、計算精度、必要な計算機資源について紹介するとともに、具体的なタンパク質をモデルとした本手法の計算例を示す。
【参考文献】[1] 平野, 千葉, 佐藤, 日本シミュレーション学会論文誌, 1, 50, (2009).
[2] F. Aquilante, L. Boman, J. Boström, H. Koch, R. Lindh, M. A. Sánchez, T. B. Pederson,

LinearScaling Techniques in Computational Chemistry and Physics, 301, 2011.
[3] Y. Okiyama, T. Nakano, K. Yamashita, Y. Mochizuki, N. Taguchi, S. Tanaka, Chem. Phys.

Lett., 490, 84 (2010).
[4] H. Harbrecht, M. Peters, R. Schneider, Appl. Num. Math., 62, 428 (2012).
[5] M. Häser, R. Ahlrichs, J. Comput. Chem., 10, 104, (1989)[6] http://www.netlib.org/scalapack/4P097OEC による S0 から S4 状態までの水分子酸化機構に関する理論的研究(三重大院工) 市野智也・三谷昌輝・吉岡泰規
Theoretical study on oxidation of water molecule from S0 to S4 states catalyzed by OEC(Mie Univ.) Tomoya Ichino・Masaki Mitani・Yasunori Yoshioka[序]

近年、太陽光エネルギーの利用の観点から光合成の仕組みに関心がもたれている。大気中の酸素分子は光エネルギー吸収で駆動する水の四電子酸化反応 (2H2O → O2 + 4H+ + 4e－)で生成される。この反応は活性中心に CaMn4Ox クラスターをもつ酸素発生複合体 (OEC)に触媒される。電子とプロトンを段階的に放出しながら S0 から S4 の 5 つの酸化状態を経由する OEC の触媒サイクル (Kok サイクル)が提案されているが、反応中間体に関するデータが不足しているため反応機構は未だ確立されていない。
酸素発生は S4→S0 遷移で起こるとされていることから、S4 状態あるいは S0 への遷移の間で O－O 結合が生成すると想定される。しかしながら、S4 状態は他の S 状態に比べ分光学的解析が困難なため、分子構造・酸化状態は明らかにされていない。さらに、S3 あるいは S4 状態で Mn(V)－oxo 結合や Mn(IV)－oxyl 結合といった高原子価 Mn－O 結合の生成が必要とされ、この活性種と Ca2+イオンに配位する水(OH－)あるいは第二配位圏にある水との間で O－O 結合が生成されると提案されているものの、基質水分子の位置が不明瞭であるためその詳細は謎である。
2011 年に報告された高解像度 X 線構造は OEC 中の水の位置を初めて明らかにし、水の酸化反応の解明に期待されている [1]。本発表では S0 から S4 状態までの分子構造と酸化状態を理論的に明らかにし、O－O 結合生成に関与する反応物を含めた S0 から S4 状態までの水の酸化反応機構を提案する。

[計算方法]

高解像度 X 線構造 (PDB ID: 3ARC)を利用し、CaMn4 クラスターを架橋する 5 つの O 原子はすべてオキソとしてモデル化し、Asp170・Glu333・Asp342・Ala344・Glu354 はギ酸イオンに、Asp61・Glu189 は酢酸イオンに、His190・His332・His337 はイミダゾールに、Tyr161 はフェノール (PhOH)に、
Gln165 はホルムアミドに、Arg357 はグアニジウムイオンに置き換えた。Ca および Mn 原子に配位する
4 つの水を含めた計 9 つの水を含め、酸化反応に関与すると想定される水二分子に対して Kok サイクルに基づき脱プロトン化を行った。電子状態計算は低スピン状態に対して broken-symmetry 法を合わせた非制限 B3LYP 法で実行した。Mn 原子には Wachters の DZ 基底を、H・C・N・O・Ca 原子には 6-31G*
基底を使用した。構造最適化は CaMn4O5 クラスター・水分子・PhOH の H 原子に対して行い、残りの原子は X 線構造の位置で固定した。すべての計算は Gaussian09 で実行した。

[結果・考察]図 1 に OEC による S0 から S4 状態までの水の酸化反応機構を模式的に示す。図中の各 S

状態の Mn 原子の酸化数に対するスピン密度の値を表 1 に示す。S0 状態の酸化状態は(Mn1, Mn2, Mn3,
Mn4) = (III, III, IV, III)として割り当てられ、もう一つの(II, III, IV, IV)型ではない。CaMn3 キューブの外側にある Mn4(III)と Ca2+イオンのそれぞれには水二分子が配位している。
S1 状態へは Mn2(III)が Mn2(IV)へ酸化され、酸化状態は(III, IV, IV, III)となる。脱プロトン化は Ca2+
イオン上の水二分子ではなく、Mn4(III)－OH2(W2)•••OH2(W6)•••OH2(W7)•••Asp61 の水素結合ネットワークを形成する W2 である。生成した OH－は Mn4(III)に配位し、その配位距離は 1.835 Å である。

図 1. OEC による S0 から S4 状態までの水の酸化反応機構. 点線は水素結合を表す.
S2 状態へは Mn4(III)が Mn4(IV)へ酸化され、酸化数は(III, IV, IV, IV)となる。Mn4原子の電荷の増大により Mn4(IV)－OH－
結合長は 1.782 Å と短くなる。
S3 状態へは Mn4(IV)に配位する OH－が脱プロトン化と一電子酸化されて oxyl ラ
ジカル種 (O•－ )となる。S2 →S3 遷移では
Mn 原子あるいは水誘導配位子の酸化を伴うとされているが、計算は Mn 原子が酸化され(IV, IV, IV, IV)型になるのではなく、脱プロトン化と共に水誘導配位子 (OH－)が酸化されると推測する。生成した Mn4(IV)－O•－結合長は 1.696 Å である。oxyl ラジカルは Mn4(IV)の不対電子と反強磁性的にカップリングしている。Mn4 原子に配位する W2 が oxyl ラジカル種まで酸化されることから、
この水分子が基質水分子の一つであると考えられる。この oxyl 原子と他方の O 原子として可能性のある Ca2+イオンに配位する水 (W3)および第二配位圏にある W6 との O 原子間距離はそれぞれ 3.212 Å と
2.304 Å であり、後者の方が著しく近い。S0→S1→S2 状態の間、O (W6)原子は O (W2)原子から 2.64－2.86
Å も離れているが、S3 状態になると接近する。
S4 状態の最適化構造を図 2 に示す。S3→S4 遷移で Mn4(IV)－O•－結合は Ca2+イオンに配位する W3 ではなく第二配位圏の W6 と反応し、O－O 結合が形成される。W6の脱プロトン化を伴った OOH 基が中性の OOH ラジカルである
こと、および酸化状態が(III, IV, IV, III)であることから、S3→S4
遷移での一電子酸化とともに Mn4(IV)から Mn4(III)へと一電子還元が起きる。O－O 結合長は酸素分子の二重結合長 (1.21 Å)より若干長い 1.318 Å である。OOH ラジカルは Mn4(III)の不対電子と反強磁性的にカップリングしている。この反応機構は、S4 状態で生成する Mn4(IV)－oxyl ラジカル種が Ca2+イオンに配位する水あ
るいは Ca2+イオンと CaMn3 キューブ外側の Mn 原子を架橋する O原子と反応する機構とは異なる [2, 3]。
図 2. S4 状態の最適化構造.
参考文献[1] Y. Umena et al., Nature, 2011, 473, 55-61. [2] E. M. Sproviero et al., J. Am. Chem. Soc., 2008, 130,3428-3442. [3] P. E. M. Siegbahn, Acc. Chem. Res., 2009, 12, 1871-1880.

4P-098
プロスタサイクリン合成酵素とプロスタグランジン H2 複合体に関する
分子動力学的研究（茨城大理）北條

貴之、渡部

良、吉村

誠慶、森

聖治

Molecular Dynamics Study of A Complex Between Prostaglandin H2
and Prostacyclin Synthase(Ibaraki Univ.) Takayuki Houjou, Ryo Watanabe, Takayoshi Yoshimura, Seiji Mori
＜緒言＞

Scheme 1 PGH2 から PGI2 への異性化反応

プロスタサイクリン合成酵素(PGIS)はプロスタグランジン H2(PGH2)からプロスタサイクリン(PGI2)への異性化反応を触媒する P450(ステロイドの合成や薬物代謝に関与している重要なヘムタンパク質)の一種であり、酸素と還元型補酵素を必要としない。この反応のメカニズ
ムは当研究室によって解明されたが、まわりのタンパク質や溶媒(水)の影響が考慮されていな
かった１）。本研究ではタンパク質や溶媒(水)の効果を考慮に入れ、PGH2 アナログと PGIS の複合体に関する分子動力学(MD)シミュレーションを行った。
＜計算方法＞

初 期 構 造 は 、 Minoxidil が 結 合 し た ヒ ト 由 来 の PGIS(PDB:3B6H(Resolution:1.62 Å)) を 使 用 2 ）。 こ の 構 造 か ら Minoxidil を 取 り 除 き 、 基 質 ア ナ ロ グ
U46619(Figure 1)の、2D NOESY 1H NMR のデータ 3）と矛盾しない配座を作成し、AMBER
プログラム(Ver.10・ff03,GAFF, TIP3P 力場)を使用し MD シミュレーション(3 ns)を行った。
得られたスナップショットのうち NOESY のデータと矛盾しない構造を選び、U46619 をPGH2 に変えて MD シミュレーション(10 ns)を行った。

Figure 1 基質アナログ U46619

Figure 2 プロスタサイクリン合成酵素(PDBID: 3B6H)

MD シミュレーション全体のエネルギーの平均値に近い構造を初期構造とし，QM 領域の全原子の位置から 14 Å 以上離れた水を取り除いて ONIOM(B3LYP:AMBER)法により，まず構
造最適化を行った。
＜結果＞

Figure 3 分子内水素結合あり

Figure 4 分子内水素結合なし(2006.72 ps)(9416 ps)

MD シミュレーション：PGH2 がヘムに配位した状
態で側鎖がフレキシブルに動く様子が見られた。
4000 ps まで PGH2 のα鎖のカルボキシル基とω鎖のヒドロキシ基との間で分子内水素結合を形成していた(Figure 3)が、シミュレーションが進むにつれて
分子内水素結合は外れていった(Figure 4)。
ONIOM 計算：MD シミュレーションの結果より，
PGH2 が分子内で相互作用していない構造を１つ選
び、ONIOM 法による構造最適化が終了した。その結果，鉄とエンドペルオキシド酸素の距離が 2.4 Åから 3.3 Å に離れた。おそらく，エンドペルオキシドのルイス塩基性が低いことが理由として考えられる。
今後残りの 2 つの最適化も行い、Fe と酸素の距離が離れた原因を詳しく考察するとともに、反応機構解

Figure 5 9416 ps の構造を初期構造とした ONIOM 法による最適化後の構造明に役立てる。
１）T. K. Yanai, S. Mori, Chem. Eur. J. 2009, 15, 4464-4473
２）C. W. Chiang, H. C. Yeh, L. H. Wang, N. L. Chan, J. Mol. Biol. 2006, 364, 266–274
３）K. H. Ruan, J. Wu, V. Cervantes, Biochemistry 2008, 47, 680-688.

4P-099

球面グリッドの回転による多中心数値積分の精度向上（東大生研） 松田潤一，平野敏行，佐藤文俊

A Improvement of Multicener Integration
with Rotated Spherical Grids(Institute of Industrial Science, The University of Tokyo)
Junichi Matsuda, Toshiyuki Hirano, Fumitoshi Sato【序】
密度汎関数（DF）計算で振動数を求めるには全エネルギーで 0.01～0.1kcal/mol の精度が必要となるが、交換相関（XC）項は、通常、数値積分で計算され、その精度が問題となる。タンパク質などの大規模生体分子の DF 計算では、XC 項の計算が計算時間の大部分を占めるため、要求精度をなるべく少ないグリッド点数で達成するのが望ましい。今回我々は、空間分割の weighting
function と球面積分グリッドとの関係に着目し、グリッド点数を増やすことなく球面グリッドの回転操作によって精度を向上させる方法を考案した。その内容と評価結果について報告する。
【空間分割と球面積分】
全空間積分である多中心積分は各原子を中心とする単中心積分へと展開され、単中心積分は、
空間分割、動径積分、球面積分、を実現するグリッドを用いて計算される。
𝐼𝐴 =𝑤𝐴 𝒓 𝐹 𝒓 𝑑𝑟 𝑑𝛺 =

𝑁𝑟

𝑤𝑟

𝑁𝛺

𝑤𝐴 𝒓 𝐹 𝒓 𝑤𝛺 − 𝛿𝐴(1)(1)において 𝐹 𝒓 = 1 としたときの球面積分は
𝑤𝐴 𝒓 𝑑𝛺 =𝑁𝛺𝑤𝐴 𝒓 𝑤𝛺 − 𝛿(2)となる。(2)は空間分割の weighting function である𝑤𝐴 の球面積分であり、誤差 𝛿 には動径積分
や 𝐹 𝒓 に起因する誤差は含まれない。数値計算における実際の 𝛿 をみるため、単純な 2 原子分子モデルを考える。球面グリッドを角度 𝜃, 𝜑で回転したときの 𝛿 の様子を図 1 に示す。ここでは標準的に使われている Becke の空間分割と Lebedev の球面グリッドを用いている。結合方向にグリッド点があるときに 𝛿 は極大値をとる。また、白い部分は 𝛿 = 0 となる角度である。つまり空間分割に対する球面グリッドの配置角は精度に影響する。次にその誤差の大きさについて考える。
𝛿 の最大値は約 3e-7 であり、このとき(2)の真値に対する精度は有効桁数 7 桁程度となる。これ
は窒素分子において同じ球面グリッドを用いて電子数を計算した時の精度に近い。つまり被積分関数を𝑤𝐴 とする球面積分において、すでに誤差は全体誤差の上限付近に達している。
【グリッドの回転方法】
本研究において考案した精度向上の方法について説明する。単中心積分における各層での球面グリッドを回転させることで(2)の 𝛿 を小さくし、これに伴って(1)の 𝛿𝐴 が減少することを狙う。
(1)の被積分関数は𝑤𝐴 𝒓 𝐹 𝒓 なので、必ずしも𝛿が小さくなるとは限らないが、𝛿𝐴 における𝑤𝐴 の影響が𝐹 𝒓 よりも十分大きければ、精度の向上が期待できる。回転角の決定には(2)の真値が必要となるが、多原子分子ではこれを求めることが難しいため、代わりに 𝑁𝛺 = 5810 の球面グリッド
での計算値を用いる。すべての球面グリッドについて、適当な初期角度からニュートン法によって回転角 𝜃, 𝜑 を求め、あらたな角度を決定する。局所解に達した場合は新たな初期角度からこれを繰り返す。
【評価結果】(1)において𝐹 𝒓 = 𝜌 𝒓 とし、分子の持つ電子数

𝐴 𝐼𝐴 を計算して、その精度を評価した。本

方式を適用した場合と、全球面グリッドを一律で回転させた場合の最悪値とを比較した。表 1 に
結果の一部を示す。表中の数値は、Acc = −log10𝛿𝐴

𝐼𝐴であり、精度の桁数を一般化したものである。アミノ酸 1 残基において常に 1 桁程度の精度向上がみられた。参考として図 2 に、Gly の計算において全球面グリッドを一律で回転させた場合の誤差分布を示す。本方式の適用によって、分布中央付近の高い精度が得られたことがわかる。
本方式は統計的に積分精度の向上を期待するものであり、実際の運用で常に効果が得られるかについて更に多くの評価を行う必要がある。原理と方法の詳細、および他のアミノ酸や規模の大きな生体分子での評価結果を当日発表する。

図 1：球面グリッドの回転に対する𝑤𝐴 𝑤𝛺 の誤差 𝛿

図 2：球面グリッドを一律回転したときの原子配置：(0,0,0)(1,0,0). 空間積分：Becke の方法.

誤差の分布(Gly)

球面積分：Lebedev グリッド 𝑁𝛺 = 194, 半径=1.3

上辺の数値は対応する Acc の値．線および 𝜎, 𝜇 の図中の文字は結合方向におけるグリッド点を示す

値は正規分布フィッティングの結果．

表 1：電子数計算での精度 Acc の比較
Gly

Ile

Val

Glu

Gln

this work

6.00

5.64

5.79

5.85

5.36

worst

4.66

4.28

4.36

4.54

4.52

4P100溶液中のバクテリオクロロフィル a の励起状態に関する理論的研究（分子研 1，学振 PD2，ワシントン大・生化 3，海外学振 4，京大院・理 5，総研大 6）

東 雅大 1,2，小杉 貴洋 3,4，林 重彦 5，斉藤 真司 1,6Theoretical Studies on Excited States of Bacteriochlorophyll a in Solutions(Institute for Molecular Science1, JSPS Research Fellow2，University of Washington3，
JSPS Postdoctoral Fellow for Research Abroad4，University of Kyoto5, SOKENDAI6）
Masahiro Higashi1,2, Takahiro Kosugi3,4, Shigehiko Hayashi5, Shinji Saito1,6

クロロフィルやバクテリオクロロフィルは光合成で重要な役割を果たす分子である。近年、単量体 1 つにつきバクテリオクロロフィル a (BChl a)を 7 つ含む光捕集アンテナの FennaMatthews-Olson タンパク質において、長時間量子コヒーレンスを伴って励起エネルギー移動が起こることが実験的に示されて注目を集めている。そこでは、異なる環境に置かれた各サ
イトの BChl a の励起エネルギーは、約 1 kcal/mol と非常に狭い範囲内に存在していると考えられている。
このような系のダイナミクスを分子レベルで取り扱うためには、BChl a の励起状態の性質を正確に取り扱う必要がある。そこで我々は、まず溶液中の BChl a の性質に注目する。
BChl a の第一励起状態への吸収スペクトルは、様々な溶媒に対してほとんど変化しないという性質があるが、何故そうなるのかよく分かっていない。本研究では、QM/MM-RWFE-SCF
法(J. Chem. Theory Comput., 2012, 8, 322)を用いて、いくつかの溶媒中で BChl a の励起エネル
ギー計算を行った。その結果、吸収スペクトルの溶媒依存性が非常に小さいのは、基底状態から第一励起状態に遷移した時の Dipole の減少による不安定化と負電荷がカルボニル基に飛ぶこと
による溶媒との水素結合の安定化が打ち消し合うことによって、溶質－溶媒間の相互作用が基底状態と励起状態でほぼ一定に保たれるためであることが分かった。また、用いる電子状態手法に結果が大きく依存することも分かった。詳細は当日議論する。